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ABSTRACT 

We  present  a  model  for  multi-wavelength  mixing  in  semiconductor  optical  amplifiers  (SOAs)  based  on  coupled-mode  equations. 
The  proposed  model  applies  to  all  kinds  of  SOA  structures,  takes  into  account  the  longitudinal  dependence  of  carrier  density 
caused  by  saturation,  it  accommodates  an  arbitrary  functional  dependencies  of  the  material  gain  and  carrier  recombination  rate 
on  the  local  value  of  carrier  density,  and  is  computationally  more  efficient  by  orders  of  magnitude  as  compared  with  the 
standard  full  model  based  on  space-time  equations.  We  apply  the  coupled-mode  equations  model  to  a  recently  demonstrated 
phase-sensitive  amplifier  based  on  an  integrated  SOA  and  prove  its  results  to  be  consistent  with  the  experimental  data.  The 
accuracy  of  the  proposed  model  is  certified  by  means  of  a  meticulous  comparison  with  the  results  obtained  by  integrating  the 
space-time  equations. 

INTRODUCTION/PROBLEM  STATEMENT 

Semiconductor  optical  amplifiers  (SOAs)  have  been  in  the  spotlight  for  many  years,  attracting  ever  growing  interest  in  multiple 
areas  of  applications.  These  include  all-optical  signal  processing  in  fiber-optic  communication  networks,  cost-effective  local  area 
transmission,  and,  more  recently,  integrated  silicon  photonics,  where  SOAs  are  the  building  blocks  for  the  implementation  of 
large-scale  integrated  photonic  circuits.  Many  of  these  applications  rely  on  the  mixing  of  the  wavelength  components  of  the 
propagating  electric  field,  and  their  theoretical  study  can  be  performed  by  numerically  integrating  the  coupled  nonlinear 
equations  describing  the  evolution  of  the  electric  field  envelope  in  the  longitudinal  direction  along  the  SOA,  and  the  temporal 
carrier  dynamics  [1,2].  Obviously,  this  approach  is  not  suitable  for  the  efficient  design  of  an  SOA,  owing  to  the  intensive 
computational  effort  that  it  involves.  The  search  for  computationally  efficient  and  analytically  tractable  models  has  yield  the 
formulation  of  what  is  sometimes  referred  to  as  a  “reduced  model”  for  the  nonlinear  SOA  response  [3],  where  the  space-time 
equations  reduce  to  a  single  ordinary  differential  equation  [3],  suitable  for  the  analytic  study  of  multi-wave  mixing  (see  e.g.  [3- 
5]).  The  formulation  of  a  reduced  model  hinges  upon  two  major  assumptions.  The  first  is  that  the  spontaneous  carrier 
recombination  rate  is  proportional  to  the  carrier  spatial  density,  and  the  second  is  that  the  material  gain  also  depends  linearly  on 
the  carrier  density.  These  assumptions  emanate  from  early  studies  of  semiconductor  lasers.  Indeed,  in  lasers  the  carrier  density 
dynamics  is  characterized  by  small  deviations  from  a  steady  state  value  which  is  set  by  the  threshold  condition  of  gain  equalling 
the  cavity  loss.  The  small  deviations  around  this  value  are  only  caused  by  amplified  spontaneous  emission  (ASE)  and  by  some 
spatial  hole  burning,  which  is  however  of  little  significance  because  in  most  structures  the  intra-cavity  optical  intensity  is  only 
moderately  inhomogeneous.  Consequently,  in  laser  structures,  gain  and  spontaneous  emission  rate  can  be  accurately 
described  by  a  linearized  expression  around  the  steady  state  carrier  density.  Early  studies  on  SOA  structures  also  used  linear 
expressions  for  gain  and  carrier  recombination,  and  in  this  case  the  linearization,  albeit  less  accurate,  found  its  ground  on  its 
simplicity  and,  more  importantly,  on  the  limited  gain  of  legacy  SOAs,  which  implied  a  limited  longitudinal  inhomogeneity  of  the 
optical  field  in  the  optical  waveguide. 

Unfortunately,  these  assumptions  do  not  reflect  the  characteristics  of  modern  SOAs,  as  is  clarified  in  what  follows.  Modern 
SOAs  may  have  linear  gain  in  excess  to  40dB,  implying  a  pronounced  longitudinal  inhomogeneity  of  the  field  intensity  and 
hence  of  gain  saturation.  This  may  cause,  in  some  cases,  that  the  gain  is  only  slightly  saturated  at  the  waveguide  input, 
whereas  it  is  almost  zero  at  the  waveguide  output,  where  saturation  is  so  high  that  the  carrier  density  approaches  its 
transparency  value.  When  this  is  the  case,  a  linear  expression  for  the  gain  is  reasonably  accurate  only  if  the  gain  does  not 
deviate  significantly  from  the  linear  expansion  around  the  transparency  carrier  density  over  a  range  of  values. 

The  nowadays  widely  accepted  forms  for  the  dependence  of  the  material  gain  on  carrier  density  do  not  meet  this  requirement, 
because  over  such  wide  range  of  carrier  density  values  the  nonlinearity  cannot  be  neglected,  especially  in  quantum-well  (QW) 
SOAs  devices  [6].  This  makes  the  use  of  linear  forms  for  the  gain  not  an  option  for  an  accurate  and  quantitative  description  of 
the  SOA  dynamics.  In  addition,  advances  in  material  fabrication  have  made  in  modern  devices  the  contribution  of  defect- 
induced  carrier  recombination,  which  is  proportional  to  the  carrier  density  N,  negligible,  with  the  consequence  that  spontaneous 
carrier  recombination  is  dominated  primarily  by  radiative  recombination,  whose  rate  is  proportional  to  NA2,  and  secondarily  by 


Auger  recombination,  whose  rate  is  proportional  to  NA3  [6].  This  reality  makes  the  linearization  of  the  spontaneous 
recombination  rate  also  a  questionable  approximation.  All  these  arguments  together  suggest  that  the  accuracy  of  models  of  the 
nonlinear  SOA  response  based  on  linearization  of  the  carrier  recombination  rate  and  gain  may  be,  in  state-of-the-art  devices, 
highly  inaccurate. 

A  natural  approach  to  the  study  of  wave  mixing  in  SOAs,  which  closely  reminds  coupled-mode  theories,  is  the  one  based  on  the 
derivation  of  evolution  equations  for  the  complex  amplitudes  of  the  field  frequency  components.  Somewhat  surprisingly,  studies 
of  wave  mixing  in  modern  SOAs  (that  is,  SOAs  characterized  by  a  nonlinear  dependence  of  the  recombination  rate  and  material 
gain  on  carrier  density)  based  on  this  approach  seem  to  be  absent  in  the  literature.  In  a  couple  of  recent  papers  [7,8],  the 
authors  assume  a  linear  gain  and  a  polynomial  recombination  rate,  as  it  would  be  appropriate  for  bulk  SOAs.  However,  they 
express  the  recombination  rate  as  R(N)  =  NAtau_c(N),  where  \tau_c(N)  =  N/R(N)  has  the  meaning  of  an  equivalent 
spontaneous  carrier  lifetime  and,  in  the  derivation  of  the  coupled-mode  equations,  they  replace  tau_c(N)  with  some  time-  and 
space-independent  value.  This  makes,  again,  the  assumed  carrier  recombination  rate  linear. 

Another  distinctive  assumption  of  all  existing  coupled-mode  approaches  to  multi-wave  mixing  in  SOAs  is  that  the  carrier  density 
modulation  induced  by  the  mixing  is  characterized  by  a  single  harmonic  component  [8].  This  is  a  reasonable  assumption  when 
a  single  frequency  component  is  dominant  over  the  others,  like  for  instance,  in  four-wave  mixing  (FWM)  experiments  where  a 
single  pump  and  a  frequency-detuned  weak  signal  are  injected  into  the  SOA.  On  the  contrary,  this  assumption  is  not  satisfied 
when  multiple  frequency  components,  detuned  by  a  few  gigahertz,  have  comparable  intensities.  This  configuration 
characterizes  for  instance  experiments  where  two  strong  pumps  are  injected  at  frequencies  -\Omega+\omega_0  and  \Omega  + 
\omega_0,  and  one  is  interested  in  the  amplification  of  a  weak  signal  injected  at  the  central  frequency  \omega_0.  In  this  case, 
the  strongest  carrier  modulation  occurs  at  the  beat  frequency  2  \Omega  between  the  two  strong  pumps,  but  the  signal 
amplification  is  mainly  affected  by  the,  possibly  weaker,  carrier  modulation  at  frequency  \Omega.  This  configuration  recently 
became  of  great  interest  because  it  describes  the  operation  of  a  relevant  class  of  SOA-based  phase  sensitive  amplifiers  (PSAs) 
[9-13], 

In  this  paper,  we  derive  coupled-mode  equations  describing  multi-wavelength  mixing  in  SOAs  characterized  by  arbitrary 
functional  dependencies  of  the  recombination  rate  and  material  gain  on  carrier  density.  These  include  both  QVV  and  bulk  SOAs. 
The  proposed  model,  which  in  what  follows  we  refer  to  as  the  “couple-mode  model,”  takes  into  account  the  frequency 
dependence  of  the  material  gain,  as  well  as  all  orders  of  the  waveguide  dispersion,  and  accommodates  input  optical  waveforms 
consisting  of  arbitrary  combinations  of  multiple  frequency  components  (We  consider  here  only  the  nonlinearity  that  comes  from 
carrier  modulation,  neglecting  ultrafast  nonlinearity  arising  from  carrier  heating,  two  photon  absorption  and  spectral  hole 
burning.  This  choice  has  been  motivated  to  keep  the  analysis  simple,  and  also  because  we  are  interested  to  cases  where 
nonlinearity  is  large  enough  to  be  used  in  all-optical  processing  applications  or  to  be  an  issue  in  applications  where  linearity  is 
sought  for.  In  these  cases,  the  frequency  detuning  does  not  exceed  a  few  tens  of  gigahertz,  and  in  this  detuning  range  the 
nonlinear  modulation  is  mostly  caused  by  carriers.  The  inclusion  of  ultrafast  processes,  however,  does  not  pose  any  conceptual 
difficulties,  and  can  be  done  along  the  lines  of  ref.  [14]  assuming  that  the  gain  depends  on  quantities  other  than  carrier  density, 
like  e.g.  the  carrier  temperature  for  carrier  heating,  or  the  energy-resolved  population  of  carriers  for  spectral  hole  burning,  and 
assuming  a  linear  decay  process  of  these  quantities  towards  their  steady  state  values.  The  effect  of  carrier  capture  and  escape 
processes  in  QW  structures  can  be  similarly  taken  into  account  by  considering  two  distinct  carrier  densities,  one  for  the 
confinement  region  and  one  for  the  QW.  Also  these  processes,  however,  become  of  relevance  for  a  pump-probe  frequency 
detuning  of  the  order  of  one  hundred  of  GHz  [15],  much  higher  than  the  range  of  values  considered  here.).  The  implementation 
of  the  model  is  illustrated  in  detail  in  the  case  of  a  QW  SOA  characterized  by  a  logarithmic  dependence  of  the  optical  gain  on 
the  carrier  density  N,  and  by  a  cubic-polynomial  carrier  recombination  rate  R(N).  The  accuracy  of  the  coupled-mode  model  is 
successfully  tested  (unlike  in  previous  related  studies)  by  means  of  a  meticulous  comparison  with  the  results  obtained  by 
integrating  the  space-time  equations  of  the  SOA  full  model.  Remarkably,  owing  to  their  inherent  simplicity,  the  coupled-mode 
equations  imply  computational  costs  by  orders  of  magnitude  smaller  than  those  required  by  the  space-time  equations,  thus 
enabling  the  efficient  characterization  of  multi-wave  mixing  in  SOA  structures,  which  would  be  otherwise  highly  impractical.  We 
then  apply  the  derived  coupled-mode  equations  to  studying  the  operation  of  a  recently  demonstrated  dual-pumped  PSA  based 
on  an  integrated  QW  SOA  [11].  We  prove  the  results  to  be  consistent  with  the  experimental  data,  and  confirm  the  excellent 
agreement  with  the  results  obtained  by  using  the  full  SOA  model. 

CONCLUSION 

To  conclude,  we  derived  a  couple-mode  model  for  multi-wave  mixing  in  SOAs  characterized  by  arbitrary  functional 
dependencies  of  the  recombination  rate  and  material  gain  on  carrier  density.  The  model  takes  into  account  the  frequency 
dependence  of  the  material  gain,  as  well  as  all  orders  of  the  waveguide  dispersion,  and  accommodates  input  fields  consisting  of 
arbitrary  combinations  of  multiple  frequency  components.  We  showed  that  the  conventional  approach  assuming  a  limited 
number  of  generated  four-wave  mixing  components  gives  inaccurate  results  when  two  waveforms  of  similar  intensities  are 
injected  into  the  SOA.  In  this  case,  our  model  gives  highly  accurate  results  if  a  sufficient  number  of  generated  components  are 
taken  into  account,  as  we  showed  by  direct  comparison  with  full  time-domain  simulations.  We  applied  the  coupled-mode  model 
to  studying  the  operation  of  a  recently  demonstrated  dual-pumped  PSA  based  on  an  integrated  QW  SOA  [11-13],  and  showed 
that  the  outcome  of  the  model  is  consistent  with  the  experimental  results. 
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Abstract — We  present  a  model  for  multiwavelength  mixing  in 
semiconductor  optical  amplifiers  (SOAs)  based  on  coupled-mode 
equations.  The  proposed  model  applies  to  all  kinds  of  SOA  struc¬ 
tures,  takes  into  account  the  longitudinal  dependence  of  carrier 
density  caused  by  saturation,  it  accommodates  an  arbitrary  func¬ 
tional  dependencies  of  the  material  gain  and  carrier  recombina¬ 
tion  rate  on  the  local  value  of  carrier  density,  and  is  computation¬ 
ally  more  efficient  by  orders  of  magnitude  as  compared  with  the 
standard  full  model  based  on  space-time  equations.  We  apply  the 
coupled-mode  equations  model  to  a  recently  demonstrated  phase- 
sensitive  amplifier  based  on  an  integrated  SOA  and  prove  its  results 
to  be  consistent  with  the  experimental  data.  The  accuracy  of  the 
proposed  model  is  certified  by  means  of  a  meticulous  comparison 
with  the  results  obtained  by  integrating  the  space-time  equations. 

Index  Terms — Nonlinear  optics,  semiconductor  optical  ampli¬ 
fiers,  wave  mixing. 

I.  Introduction 

EMICONDUCTOR  optical  amplifiers  (SOAs)  have  been 
in  the  spotlight  for  many  years,  attracting  ever  growing  in¬ 
terest  in  multiple  areas  of  applications.  These  include  all-optical 
signal  processing  in  fiber-optic  communication  networks,  cost- 
effective  local  area  transmission,  and,  more  recently,  integrated 
silicon  photonics,  where  SOAs  are  the  building  blocks  for 
the  implementation  of  large-scale  integrated  photonic  circuits. 
Many  of  these  applications  rely  on  the  mixing  of  the  wave¬ 
length  components  of  the  propagating  electric  field,  and  their 
theoretical  study  can  be  performed  by  numerically  integrating 
the  coupled  nonlinear  equations  describing  the  evolution  of  the 
electric  field  envelope  in  the  longitudinal  direction  along  the 
SOA,  and  the  temporal  carrier  dynamics  [1],  [2],  Obviously, 
this  approach  is  not  suitable  for  the  efficient  design  of  an  SOA, 
owing  to  the  intensive  computational  effort  that  it  involves.  The 
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search  for  computationally  efficient  and  analytically  tractable 
models  has  yield  the  formulation  of  what  is  sometimes  referred 
to  as  a  reduced  model  for  the  nonlinear  SOA  response  [3],  where 
the  space-time  equations  reduce  to  a  single  ordinary  differen¬ 
tial  equation  [3],  suitable  for  the  analytic  study  of  multi-wave 
mixing  (see,  e.g.,  [3]— [5]).  The  formulation  of  a  reduced  model 
hinges  upon  two  major  assumptions.  The  first  is  that  the  spon¬ 
taneous  carrier  recombination  rate  is  proportional  to  the  carrier 
spatial  density,  and  the  second  is  that  the  material  gain  also  de¬ 
pends  linearly  on  the  carrier  density.  These  assumptions  emanate 
from  early  studies  of  semiconductor  lasers.  Indeed,  in  lasers  the 
carrier  density  dynamics  is  characterized  by  small  deviations 
from  a  steady  state  value  which  is  set  by  the  threshold  condition 
of  gain  equalling  the  cavity  loss.  The  small  deviations  around 
this  value  are  only  caused  by  amplified  spontaneous  emission 
(ASE)  and  by  some  spatial  hole  burning,  which  is  however  of 
little  significance  because  in  most  structures  the  intra-cavity 
optical  intensity  is  only  moderately  inhomogeneous.  Conse¬ 
quently,  in  laser  structures,  gain  and  spontaneous  emission  rate 
can  be  accurately  described  by  a  linearized  expression  around 
the  steady  state  carrier  density.  Early  studies  on  SOA  structures 
also  used  linear  expressions  for  gain  and  carrier  recombination, 
and  in  this  case  the  linearization,  albeit  less  accurate,  found 
its  ground  on  its  simplicity  and,  more  importantly,  on  the  lim¬ 
ited  gain  of  legacy  SOAs,  which  implied  a  limited  longitudinal 
inhomogeneity  of  the  optical  field  in  the  optical  waveguide. 

Unfortunately,  these  assumptions  do  not  reflect  the  character¬ 
istics  of  modern  SOAs,  as  is  clarified  in  what  follows.  Modern 
SOAs  may  have  linear  gain  in  excess  to  40  dB,  implying  a  pro¬ 
nounced  longitudinal  inhomogeneity  of  the  field  intensity  and 
hence  of  gain  saturation.  This  may  cause,  in  some  cases,  that  the 
gain  is  only  slightly  saturated  at  the  waveguide  input,  whereas 
it  is  almost  zero  at  the  waveguide  output,  where  saturation  is  so 
high  that  the  carrier  density  approaches  its  transparency  value. 
When  this  is  the  case,  a  linear  expression  for  the  gain  is  rea¬ 
sonably  accurate  only  if  the  gain  does  not  deviate  significantly 
from  the  linear  expansion  around  the  transparency  carrier  den¬ 
sity  over  a  range  of  values.  The  nowadays  widely  accepted  forms 
for  the  dependence  of  the  material  gain  on  carrier  density  do  not 
meet  this  requirement,  because  over  such  wide  range  of  carrier 
density  values  the  nonlinearity  cannot  be  neglected,  especially 
in  quantum- well  (QW)  SOAs  devices  [6].  This  makes  the  use 
of  linear  forms  for  the  gain  not  an  option  for  an  accurate  and 
quantitative  description  of  the  SOA  dynamics.  In  addition,  ad¬ 
vances  in  material  fabrication  have  made  in  modern  devices  the 
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contribution  of  defect-induced  carrier  recombination,  which  is 
proportional  to  the  carrier  density  N,  negligible,  with  the  con¬ 
sequence  that  spontaneous  carrier  recombination  is  dominated 
primarily  by  radiative  recombination,  whose  rate  is  proportional 
to  N 2,  and  secondarily  by  Auger  recombination,  whose  rate  is 
proportional  to  N3  [6],  This  reality  makes  the  linearization  of 
the  spontaneous  recombination  rate  also  a  questionable  approx¬ 
imation.  All  these  arguments  together  suggest  that  the  accuracy 
of  models  of  the  nonlinear  SOA  response  based  on  linearization 
of  the  carrier  recombination  rate  and  gain  may  be,  in  state-of- 
the-art  devices,  highly  inaccurate. 

A  natural  approach  to  the  study  of  wave  mixing  in  SOAs, 
which  closely  reminds  coupled-mode  theories,  is  the  one  based 
on  the  derivation  of  evolution  equations  for  the  complex  am¬ 
plitudes  of  the  field  frequency  components.  Somewhat  surpris¬ 
ingly,  studies  of  wave  mixing  in  modern  SOAs  (that  is,  SOAs 
characterized  by  a  nonlinear  dependence  of  the  recombination 
rate  and  material  gain  on  carrier  density)  based  on  this  ap¬ 
proach  seem  to  be  absent  in  the  literature.  In  a  couple  of  re¬ 
cent  papers  [7],  [8],  the  authors  assume  a  linear  gain  and  a 
polynomial  recombination  rate,  as  it  would  be  appropriate  for 
bulk  SOAs.  However,  they  express  the  recombination  rate  as 
R(N)  =  N/tc(N),  where  tc(N)  =  N/R(N)  has  the  meaning 
of  an  equivalent  spontaneous  carrier  lifetime  and,  in  the  deriva¬ 
tion  of  the  coupled-mode  equations,  they  replace  tc(N )  with 
some  time-  and  space-independent  value.  This  makes,  again, 
the  assumed  carrier  recombination  rate  linear. 

Another  distinctive  assumption  of  all  existing  coupled-mode 
approaches  to  multi-wave  mixing  in  SOAs  is  that  the  carrier 
density  modulation  induced  by  the  mixing  is  characterized  by  a 
single  harmonic  component  [8],  This  is  a  reasonable  assumption 
when  a  single  frequency  component  is  dominant  over  the  oth¬ 
ers,  like  for  instance,  in  four-wave  mixing  (FWM)  experiments 
where  a  single  pump  and  a  frequency-detuned  weak  signal  are 
injected  into  the  SOA.  On  the  contrary,  this  assumption  is  not 
satisfied  when  multiple  frequency  components,  detuned  by  a 
few  gigahertz,  have  comparable  intensities.  This  configuration 
characterizes  for  instance  experiments  where  two  strong  pumps 
are  injected  at  frequencies  —ft  +  ujq  and  f!  +  u)q,  and  one  is 
interested  in  the  amplification  of  a  weak  signal  injected  at  the 
central  frequency  loq.  In  this  case,  the  strongest  carrier  modu¬ 
lation  occurs  at  the  beat  frequency  20  between  the  two  strong 
pumps,  but  the  signal  amplification  is  mainly  affected  by  the, 
possibly  weaker,  carrier  modulation  at  frequency  O.  This  con¬ 
figuration  recently  became  of  great  interest  because  it  describes 
the  operation  of  a  relevant  class  of  SOA-based  phase  sensitive 
amplifiers  (PSAs)  [9]-[13], 

In  this  paper,  we  derive  coupled-mode  equations  describing 
multi-wavelength  mixing  in  SOAs  characterized  by  arbitrary 
functional  dependencies  of  the  recombination  rate  and  material 
gain  on  carrier  density.  These  include  both  QW  and  bulk  SOAs. 
The  proposed  model,  which  in  what  follows  we  refer  to  as  the 
couple-mode  model,  takes  into  account  the  frequency  depen¬ 
dence  of  the  material  gain,  as  well  as  all  orders  of  the  wave¬ 
guide  dispersion,  and  accommodates  input  optical  waveforms 
consisting  of  arbitrary  combinations  of  multiple  frequency 


components.1  The  implementation  of  the  model  is  illustrated 
in  detail  in  the  case  of  a  QW  SOA  characterized  by  a  log¬ 
arithmic  dependence  of  the  optical  gain  on  the  carrier  den¬ 
sity  N,  and  by  a  cubic-polynomial  carrier  recombination  rate 
R(N).  The  accuracy  of  the  coupled-mode  model  is  success¬ 
fully  tested  (unlike  in  previous  related  studies)  by  means  of  a 
meticulous  comparison  with  the  results  obtained  by  integrating 
the  space-time  equations  of  the  SOA  full  model.  Remarkably, 
owing  to  their  inherent  simplicity,  the  coupled-mode  equations 
imply  computational  costs  by  orders  of  magnitude  smaller  that 
those  required  by  the  space-time  equations,  thus  enabling  the 
efficient  characterization  of  multi-wave  mixing  in  SOA  struc¬ 
tures,  which  would  be  otherwise  highly  impractical.  We  then 
apply  the  derived  coupled-mode  equations  to  studying  the  oper¬ 
ation  of  a  recently  demonstrated  dual-pumped  PSA  based  on  an 
integrated  QW  SOA  [11].  We  prove  the  results  to  be  consistent 
with  the  experimental  data,  and  confirm  the  excellent  agreement 
with  the  results  obtained  by  using  the  full  SOA  model. 


II.  Coupled-Mode  Equations  for  Multi- Wavelength 
Propagation  in  SOAs 


We  denote  by  E(z,t)  the  slowly-varying  complex  envelope 
of  the  electric  field  propagating  in  the  SOA  in  the  temporal 
reference  frame  that  accommodates  the  field  group  velocity  vg , 
corresponding  to  the  real  field 


£{z,  t)  =  Re 


E 


e 


■i[u0t-l3(u)o)z] 


(1) 


with  u) o  being  the  optical  frequency.  The  field  envelope  E  is 
normalized  so  that  that  \E\2  is  the  optical  power  flowing  through 
the  transverse  waveguide  section.  It  is  related  to  the  photon  flux 
P  in  photons  per  unit  time  and  area  through  the  relation 


\E\2  =  feVoSmod^, 


(2) 


where  S'Tnod  =  S/T  is  the  modal  area  of  the  waveguide,  with  S 
denoting  the  effective  SOA  area  and  T  the  optical  confinement 
factor.  The  evolution  of  E(z,t)  along  the  SOA  is  governed  by 
the  familiar  equation 

dE  1 

—  =  -[(1- ia)Fg  -  ctint}E +  ij3E +  rsp,  (3) 


1  We  consider  here  only  the  nonlinearity  that  comes  from  carrier  modulation, 
neglecting  ultrafast  nonlinearity  arising  from  carrier  heating,  two  photon  ab¬ 
sorption  and  spectral  hole  burning.  This  choice  has  been  motivated  to  keep  the 
analysis  simple,  and  also  because  we  are  interested  to  cases  where  nonlinearity 
is  large  enough  to  be  used  in  all-optical  processing  applications  or  to  be  an 
issue  in  applications  where  linearity  is  sought  for.  In  these  cases,  the  frequency 
detuning  does  not  exceed  a  few  tens  of  gigahertz,  and  in  this  detuning  range  the 
nonlinear  modulation  is  mostly  caused  by  carriers.  The  inclusion  of  ultrafast 
processes,  however,  does  not  pose  any  conceptual  difficulties,  and  can  be  done 
along  the  lines  of  ref.  [14]  assuming  that  the  gain  depends  on  quantities  other 
than  carrier  density,  like,  e.g.,  the  carrier  temperature  for  carrier  heating,  or  the 
energy-resolved  population  of  carriers  for  spectral  hole  burning,  and  assuming 
a  linear  decay  process  of  these  quantities  towards  their  steady  state  values.  The 
effect  of  carrier  capture  and  escape  processes  in  QW  structures  can  be  similarly 
taken  into  account  by  considering  two  distinct  carrier  densities,  one  for  the  con¬ 
finement  region  and  one  for  the  QW.  Also  these  processes,  however,  become  of 
relevance  for  a  pump-probe  frequency  detuning  of  the  order  of  100  GHz  [15], 
much  higher  than  the  range  of  values  considered  here. 
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where  a  is  the  Henry  factor,  «int  is  the  SOA  internal  loss  co¬ 
efficient,  and  rsp  is  the  spontaneous  emission  noise  term.  By  g 
and  j3  we  denote  the  material  gain  operator  and  the  wavenum¬ 
ber  operator.  The  operator  formalism  allows  us  to  conveniently 
accommodate  the  frequency  dependence  of  the  gain  as  well  as 
the  waveguide  dispersion  to  any  order.  Within  this  formalism 
the  two  operators  can  be  expressed  as 


-  _  y-  1  dng(N,  wo)  (,d\ 
9  ^  n\  dufi  Ydt) 

m.  =  f)  u  x  7 


(4) 


;  y'  1  d'llM 

P  ^  n\  d 

m  =  2  0 


(5) 


where  g(N ,  ui)  is  the  gain  coefficient  expressed  as  a  function  of 
the  carrier  density  N  and  the  optical  frequency  lu,  and  (3(u>)  is 
the  frequency-dependent  field  propagation  constant.  The  expres¬ 
sions  for  g  and  f3  in  Eqs.  (4)  and  (5)  are  obtained  by  expanding 
g(N:uj)  and  /3(w)  around  the  carrier  frequency  The  fact 
that  the  sum  in  Eq.  (5)  starts  from  n  =  2  is  consistent  with  the 
definition  of  E(z ,  t)  in  Eq.  (1),  which  already  accounts  for  the 
effect  of  /3(u>o)  and  d/3(wo)/da20  =  l/us. 

The  spontaneous  emission  noise  term  rsp  is  modeled  as  a 
zero-mean,  complex  phase  independent  random  process.  It  de¬ 
pends  explicitly  on  the  carrier  density,  besides  time  and  space, 
i.e.,  rs p  =  rSp(N7  f;  z).  Its  correlation  function  is  [16] 

E  KP z)rsp(N,t'-,  z')\  =  hw0Rsp(N,t  -  t')6(z  -  z'), 

(6) 

where  by  the  symbol  E  we  denote  ensemble  averaging.  Here  the 
term  S(z  —  z')  accounts  for  the  fact  that  different  longitudinal 
waveguide  sections  provide  statistically  independent  contribu¬ 
tions  to  the  noise  term,  and  [6] 


Rsp (N,  t'  —  t)  =  J  e-i(— o)(t'-t)nsp(iV,a;)r5(iV!u;)^|# 

7F(7) 

is  the  spontaneous  emission  rate  into  the  waveguide  mode  and  in 
the  field  propagation  direction,  with  nsp  denoting  the  population 
inversion  factor  [6].  Spontaneous  emission  is  a  small  perturba¬ 
tion  of  the  propagating  field,  so  that  we  may  safely  replace  N 
with  its  temporal  average,  thus  neglecting  the  effect  of  its  small 
fluctuations  around  this  value.  Within  this  approximation  the 
process  of  spontaneous  emission  can  be  modeled  as  a  stationary 
process  in  time. 

The  equation  for  the  carrier  density  is 

~KT  =  Rj  ~  -Rrad  ~  Rm ,  (8) 

at 

where  the  meaning  of  each  of  the  terms  at  the  right-hand  side 
of  the  equation  is  discussed  in  what  follows.  The  term 


Rj 


JwaL 

V 


J 

ed’ 


(9) 


is  the  carrier  injection  rate  into  the  active  volume  V  =  SL  = 
wadL,  where  wa  is  the  active  region  width,  L  is  the  active 
region  length,  and  d  is  the  active  region  thickness.  The  term 
Rvad  =  Rin  +  -Rout  is  the  radiative  recombination  rate  related  to 


processes  in  which  the  recombination  of  one  carrier  is  associated 
to  the  generation  of  one  photon.  The  term  Rul  refers  to  processes 
in  which  emission  occurs  into  the  guided  mode.  By  definition, 
this  implies  that  i?jn  (TV)  is  related  to  the  flux  of  photons  flowing 
in  the  waveguide  P  through  the  balance  relation 

RinSdz  =  Sj„od  [P(z  +  dz)  -  P(z)} ,  (10) 


which  yields 


R'u 


1  dP 

F  dT 


(ID 


The  processes  accounted  for  by  RlVL  includes  stimulated  emis¬ 
sion  and  spontaneous  emission  within  the  waveguide  mode 
(which  is  only  a  fraction  of  the  overall  spontaneous  emission). 
The  term  Rout  (N)  is  the  rate  of  recombinations  accompanied  by 
spontaneous  emission  of  photons  outside  the  waveguide  mode 
and  it  can  be  expressed  as  i?out  =  BN2  —  /?s]l  ill ,  where  the 
term  BN2  is  known  to  be  an  excellent  approximation  of  the 
total  rate  of  recombinations  associated  with  spontaneous  emis¬ 
sion  (inside  and  outside  the  waveguide  mode)  [6],  and  /isp  ill 
accounts  for  the  rate  of  recombinations  that  produce  sponta¬ 
neous  emission  in  the  waveguide  mode,  namely  it  accounts  for 
the  contribution  to  the  carrier  recombination  rate  of  the  noise 
term  rsp  that  appears  in  Eq.  (3).  Finally,  the  term  RnI  is  the 
recombination  rate  of  non-radiative  processes,  which  we  ex¬ 
press  as  Rnl  =  AN  +  CN3 ,  where  the  linear  contribution  AN 
is  mostly  due  to  defect-induced  recombination,  and  the  cubic 
contribution  CN3  to  Auger  recombination.2  By  combining  the 
various  mechanisms,  Eq.  (8)  becomes 


dN 

Idt 


—R{N) 


ed  +  Rsp’in 


1  d\E\2 

HujqS  dz 


(12) 


where  by  R(N)  we  denote  the  familiar  recombination  rate3 


R(N)  =  AN  +  BN2  +  CN3.  (13) 


By  expanding  the  derivative  d\E\2/dz  and  using  Eq.  (3), 
Eq.  (12)  assumes  the  form 

BN  J  1 

_  =  -R(JV)+--_Re[£*(1-»)rSi=],  (14) 


where  we  used  the  fact  that  (3  is  a  Hermitian  operator,  and  hence 
it  does  not  contribute  to  d\E\2 /dz.  The  last  term  at  the  right- 
hand  side  of  Eq.  ( 14)  reduces  to  the  familiar  form  Tg  |  E\ 2  / hto 0  S 
if  the  gain  coefficient  is  assumed  to  be  frequency  independent. 
Note  that  in  Eq.  (14),  the  term  Rsp,m  disappears  because  it 


2We  note  that,  while  the  resulting  cubic  polynomial  expression  AN  + 
BN 2  +  CAT3  has  been  shown  to  fit  very  well  the  experimental  data  in  most 
cases  [6],  the  one-to-one  correspondence  between  the  three  terms  of  the  poly¬ 
nomial  and  the  three  recombination  mechanisms  is  not  always  as  definite  as  is 
illustrated  in  the  main  text.  For  instance,  in  the  case  of  non-parabolic  bands  (the 
normal  case),  radiative  recombination  is  also  non-parabolic  and  is  best  modeled 
with  a  bit  of  linear  component;  carrier  leakage  (due  to  finite  QW  barriers)  has  an 
exponential  dependence  and  requires  a  polynomial  fit,  affecting  the  numerical 
values  of  A,  B,  and  C. 

3 This  expression  of  R(N)  is  widely  established  and  is  given  here  for  consis¬ 
tency  with  previous  studies.  We  stress,  however,  that  the  analysis  that  follows 
does  not  make  use  of  it  explicitly,  and  rather  applies  to  arbitrary  expressions  of 

R(N). 
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cancels  with  an  opposite  term  that — by  definition — comes  from 

d\E\2 /dz. 

We  note  that  Eqs.  (3)  and  (14)  can  be  generalized  so  as  to 
include  the  field  polarization  in  the  analysis.  While  this  task 
is  rather  straightforward  and  does  not  involve  any  conceptual 
challenge,  we  intentionally  ignore  polarization-related  issues  in 
order  to  keep  the  focus  on  the  main  objective  of  this  work,  which 
is  the  study  of  multi-wavelength  propagation. 

We  express  the  multi-wavelength  electric  field  and  the  carrier 
density  as  follows, 

E(z,t)  =  (15) 

k 

N(z,t)  =  N0(z)+^ANk(z)e-iknt,  (16) 

k 

where  the  coefficients  ANk  must  satisfy  the  equality 
A N-k(z)  =  A N£(z)  for  N{z,  t)  to  be  real.  The  term  N0{z)  + 
ANo(z)  is  the  z-dependent  time-independent  value  of  the  car¬ 
rier  density  that  characterizes  the  system  when  it  achieves  its 
stationary  state,  and  Nq  (z)  is  defined  as  the  solution  of 

J-  =  R(N0)+-]—y2g(N0,uJk)\Ek\2.  (17) 

ed  fUdn  D 

k 


and  [8].  The  difference  between  these  two  quantities  is  approxi¬ 
mately  a  factor  of  2  when  the  radiative  bimolecular  recombina¬ 
tion  BN 2  is  the  dominant  contribution  to  R(N),  or  3  when  the 
Auger  recombination  CN 3  is  dominant.  By  inserting  Eqs.  (18) 
and  (19)  into  Eq.  (3)  and  Eq.  (14)  we  obtain 


dE  1 

—  =  -  [(1  -  ia)T(go  +  A NgN)  -  aint]  E  +  i(3E  +  r,  (22) 
oz  z 

and 


dAN 

dt 


AN 


R(Wo)  -  2 


t(N0) 

Re  [ E*{l-ia)Tg0E } 


-AN 


HijJqS 

Re  [E*{  1  —  ia)Tg^E] 
HloqS 


(23) 


where  the  operators  go  and  are  defined  as  in  Eq.  (4),  pro¬ 
vided  that  g(N,ui)  is  replaced  with  g{No,u>)  and  5,v (No,uj), 
respectively. 

The  evolution  equation  for  the  electric  field  coefficient  Ek  is 
obtained  by  inserting  the  expression  of  the  field  (15)  in  Eq.  (22) 
and  by  equating  the  coefficient  of  the  term  exp(—ikfli)  at  the 
two  sides  of  the  resulting  equation.  As  a  result,  one  finds 


For  values  of  the  frequency  spacing  f l /2tt  that  exceed  the  SOA 
modulation  bandwidth,  the  temporal  fluctuations  of  N(z,  t) 
around  its  stationary  value  are  filtered  by  the  carrier  dynamics 
and  hence  they  can  be  treated  within  a  perturbation  approach. 
A  consequence  of  this  situation  is  that  the  deviation  ANo(z) 
of  the  stationary  carrier  density  value  from  Nq(z)  is  also  a  per¬ 
turbation,  and  is  small  compared  to  N(){z).  In  this  framework, 
all  carrier-density  dependent  quantities  that  appear  in  Eqs.  (3) 
and  (14)  can  thus  be  expanded  to  first  order  with  respect  to 
AN  =  N  —  Nq,  namely 

AN 

R(N)=  R(N0)  +  — — ,  (18) 

r{N0) 

g(N,u)=  g(N0,u)  +  gN(N0,u)AN,  (19) 


where  by  the  subscript  N  we  denote  differentiation  with  respect 
to  N.  The  quantity 


t{N0)  =  Rn(N0 ) 


1  _ 

d  R{N) 

dN 

N=N0 

is  the  spontaneous  carrier  lifetime,  and 


gN  {N0 ,  ui) 


dg(N,u) 


dN 


N  —  No 


(20) 


(21) 


is  the  differential  gain.  We  stress  that  these  are  z-dependent 
quantities,  owing  to  the  fact  that  Nq  =  Nq{z),  and  hence  their 
values  evolve  along  the  SOA.  We  also  notice  that  the  effective 
carrier  lifetime  governing  the  dynamics  of  carrier  modulation 
around  the  steady  state  value  is  the  differential  carrier  lifetime 
t(jVo)  given  in  Eq.  (20)  and  also  introduced  in  [17],  and  not 
the  total  carrier  lifetime  tc{Nq)  =  Nq/R(Nq )  used  in  Refs.  [7] 


d  Ek 
dz 


-(1  -  ia)Tg{No,uJk)  -  ctmt  +  if)(u>k) 


Ek 


+  -(1  -  ia)  ^2  ANk_nTgN{No,Lun)En  +  rk,( 24) 


where  we  used  goE  =  di^o ,uk)Ek  exp(-ikflt)  and  g n 
E  =  gi<r{No, wk)Ek  exp(— ikflt),  with  utk  =  +  kEl. 

The  noise  term  rk  is  defined  by 


rk  (N;  z) 


d te'kMrBp(N,t,  z ) 


(25) 


has  zero  mean  ( rk(N ;  z))  =  0,  and  its  variance  follows  from 
(■ r*k{N ;  z)rh(N;  z'))  =  5{z  -  z')hio 0 

x  J  dt  j  df  exp[i£l{kt'  —  ht)\Rep{N,  f  —  t).  (26) 

Using  the  stationarity  of  Rsp,  we  may  express  the  above  as 

{r*k {N-  z)rh  {N-  z'))  =  6k,hS{z  -  z')huj0 

xnsp(JV,  wo  +  k£i)Tg(N,  uto  +  kfl).  (27) 

The  terms  rk{N\  z'),  k  =  0,  ±1,  ±2,  . . .  are  therefore  a  set  of 
phase-independent,  spatially-uncorrelated  noise  terms,  which 
can  be  modeled  as  differentials  of  independent  Wiener  pro¬ 
cesses.  At  this  point  we  can  recast  Eq.  (24)  in  the  following 
compact  form 


dE 
d z 


—  (1  —  io;)r(G  +  H)  —  aint  I  +  *b 


E- 


(28) 


where  E  and  r  are  column  vectors  constructed  by  stacking  the 
electric  field  coefficients  Ek  and  the  noise  projections  rk  one  on 
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top  of  another,  respectively,  with  Eo  and  ?’o  occupying  the  cen¬ 
tral  position,  namely  E  =  [...,  Ei ,  E\,  Eq,  E-i,  £?_ 2,  ■  •  .]*, 
and  the  same  for  r  (the  superscript  t  stands  for  “trans¬ 
posed”).  The  vector  E  and  r  are  of  course  infinite-dimensional, 
and  so  are  the  square  matrices  G,  H  and  b.  Consistently 
with  the  definition  of  E,  we  use  positive  and  negative  in¬ 
dices  to  identify  the  elements  of  these  matrices,  with  the 
(0,0)  element  occupying  the  central  position.  In  particu¬ 
lar,  G  and  b  are  diagonal  matrices  whose  (k,  k )  elements 
are  equal  to  Gk,k  =  g(N0,uk )  and  bktk  =  j3{uk)  -  (3(u0)  - 
kQd/3(u>o) /dojQ,  respectively,  whereas  the  (k,  n)  element  of  H 
isHktTl  =  ANk~ng]y(No:  un).  By  I  we  denote  the  identity  ma¬ 
trix  (regardless  of  its  dimensions). 

We  now  proceed  to  the  extraction  of  the  carrier  density  coeffi¬ 
cients  A Nk  by  equating  the  terms  proportional  to  exp(— ikilt) 
at  the  two  sides  of  Eq.  (23),  when  the  expression  of  AN  in 
Eq.  (16)  is  inserted  in  it.  After  some  straightforward  algebra, 
involving  the  use  of  Eq.  (17),  one  obtains 

(l-ikTn)ANk  =  -Y,^NhPk,h+Afk,  (29) 

h 

where 


A4  =  —t(Nq)R(No){1  —  <5fc,  0) 
'(1  -  ia)En+kE* 


E 


Pk,h 


=  E 


-fstim  (N0  ,  ^ n  +  k ) 
(1  -  ia)En+k hE* 


-^sat  (A^o  ,  ^n+k— h) 


(1  +  ia)En+kE* 


-fstim  (N0,  U)n  ) 

(1  +  ia)En+k hE\ 


,(30) 


-fsat(A^),  ^n) 


(31) 


The  quantity 


-Fsat(Ao,  w) 


HujqS 

T(N0)rgN(N0,io) 


(32) 


is  the  familiar  saturation  power,  although  its  definition  accounts 
for  the  frequency  dependence  of  the  gain  coefficient  explicitly, 
and 


Pstim{N0,Uj) 


R(N0) 


hiooS 
r  g(N0,ui) 


(33) 


is  the  power  value  above  which  carrier  depletion  is  dominated 
by  stimulated  emission.  We  hence  refer  to  Pstim  as  to  stimulated 
power.  Equation  (29)  can  be  conveniently  recast  in  the  following 
compact  form 

(I  -  irClk  +  p)AN  =  J\f,  (34) 


where  the  (k,  h)  element  of  the  matrix  p  is  equal  topkk,  and  k  is 
a  diagonal  matrix  with  diagonal  elements  nk  k  =  k.  The  column 
vectors  AN  and  N  are  constructed  (like  the  field  vector  E)  by 
stacking  the  coefficients  A  Nk  and  Afk  one  on  top  of  another,  re¬ 
spectively,  namely,  AN  =  [. . . ,  AN\,  ANq,  AiV_i,  . .  .]*  and 
Af=  [. . . ,  A/j ,  0,  AC-i ,  . . .]; .  The  coefficients  Nq  and  A Nk  are 
hence  obtained  for  a  given  electric  field  state  by  solving  Eqs.  (17) 
and  (34).  These  are  the  most  general  coupled-mode  equations 
accounting  for  any  functional  dependence  of  the  recombination 


rate  and  material  optical  gain  on  carrier  density,  as  well  as  for 
the  frequency  dependence  of  the  gain  and  waveguide  dispersion. 

III.  Implementation  of  the  Coupled-Mode  Equation 
Model  in  Realistic  SOA  Structures 

As  is  customarily  done  in  most  studies  of  practical  relevance, 
where  the  waveguide  dispersion  and  the  frequency  dependence 
of  the  gain  coefficient  have  been  shown  to  play  a  minor  role,  in 
this  section  we  neglect  chromatic  dispersion,  as  well  as  higher- 
order  dispersion,  and  assume  frequency-independent  gain.  With 
this  simplification  the  matrices  G,  H,  and  p  become  frequency- 
independent  and  assume  a  very  convenient  form,  as  is  shown 
in  what  follows.  We  also  neglect  the  presence  of  spontaneous 
emission  noise  terms,  whose  implications  on  the  SOA  perfor¬ 
mance,  chiefly  on  the  SOA  noise  figure,  will  be  the  subject  of 
future  work. 

The  multi-wavelength  propagation  model  introduced  in  the 
previous  section  involves  an  infinite  number  of  coefficients  Ek 
and  A Nk,  a  situation  that  is  obviously  incompatible  with  its 
implementation  in  any  numerical  platform.  However,  as  will 
be  shown  in  the  next  section,  high-order  coefficients  (namely 
Ek  and  ANk  coefficients  with  large  values  of  k  \ )  provide  a 
negligible  contribution  to  the  solution  of  Eqs.  (17),  (28),  and 
(34),  and  hence  they  can  be  omitted  by  truncating  the  vectors 
E  and  AN.  The  truncation  of  E  and  AN  requires  of  course 
that  all  matrices  involved  in  Eqs.  (28)  and  (34)  be  also  truncated 
accordingly.  In  what  follows  we  provide  explicit  expressions  for 
those  matrices  and  discuss  the  procedure  that  allows  the  efficient 
computation  of  E  and  AN. 

The  truncation  procedure  of  the  infinite  set  of  equations  (28) 
can  be  performed  in  a  number  of  ways.  One  possible  approach 
is  assuming  that  Ek(z)  =  0  for  |fc|  >  M.  Here  M  is  an  integer 
number  that  can  be  determined  self  consistently  by  checking 
that  the  integration  of  the  equations  for  M  — >  M  +  1  yields  in¬ 
distinguishable  results.  This  assumption  implies  A Nk(z)  =  0 
for  |fc|  >  2 M,  owing  to  the  absence  of  beat  terms  at  frequency 
offsets  larger  than  2Mfi.  A  simpler  yet  equally  accurate  ap¬ 
proach  is  to  assume  that  the  carrier  density  coefficients  ANk  (z) 
are  also  zero  at  frequency  offsets  greater  than  Mil.  Here  we 
adopt  the  latter  approach,  within  which  Eqs.  (15)  and  (16) 
specialize  to 

M 

E(z,t)  =  ]T  Ek(z)e~imt,  (35) 

k  —  —M 

M 

N(z,t)=  N0(z)+  A Nk(z)e~lkSlt.  (36) 

k=-M 

Accordingly,  the  field  vector  E  and  carrier  density  modula¬ 
tion  vector  AN,  consist  of  (2 M  +  1)  components.  Matrices 
G  and  H  in  Eq.  (28)  become  (2 M  +  1)  x  (2 M  +  1)  matri¬ 
ces.  In  particular,  owing  to  the  assumption  of  frequency-flat 
gain,  one  can  readily  verify  the  equalities  G  =  g(No)l,  and 
H  =  </iv(iVo)T(AiV),  where  by  T2M+1  (AA^)  we  denote  a 
Hermitian- symmetric  Toeplitz  matrix  [18].  Below  we  give  the 
expression  of  T2M+i(ANk)  in  the  case  M  =  2  for  illustration 
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purposes, 


T5(AiVfc) 


AA0 

A  Nj_ 

A  N2 

0 

0 

AN£ 

A  N0 

ANi 

an2 

0 

AN 2* 

ANi 

A  No 

ANi 

an2 

0 

AN 

AN* 

AN0 

ANi 

0 

0 

AN% 

A  AT* 

AN0 

The  neglect  of  the  waveguide  dispersion  yields  b  =  0,  and  hence 
Eq.  (28)  simplifies  to 


d  E 
cb 


(1  -  ia)g(N0)  -  aint 
2 


I  +  T2M+i(A7Vit) 


E,  (38) 


where  Nq  is  the  solution  of 


J 

ed 


R(N0 ) 


1  + 


\E? 

T-Pstim(iVo) 


(39) 


with 

Pst^(N0)  =  R(N0) 


(40) 


The  expression  for  the  carrier  density  modulation  vector  AN 
simplifies  to 


A  N=- 


tR{N0 ) 

-Pstim(Ao) 


I  -  rOk  + 


T2M  +  l(Cfc) 

PsM) 


-1 


c, 


(41) 


TABLE  I 

SOA  Parameters 


Description 

Value 

Units 

Lineal'  recombination  coefficient  A 

106 

s_1 

Bimolecular  recombination  coefficient  B 

0.3  x  10“10 

cm3  /  s 

Auger  coefficient  C 

3.3  x  10~29 

cm6  /s 

Optical  confinement  factor  T 

9.7% 

Linewidth  enhancement  factor  a. 

5 

Optical  wavelength  Ao 

1561 

nm 

Group  velocity  vg 

8.33  x  109 

cm/s 

Active  region  width  wa 

2  x  10~4 

cm 

Active  region  tickness  d 

65  x 10~7 

cm 

Active  region  length  L 

0.1 

cm 

Gain  coefficient  go 

1800 

cm-1 

Transparency  carrier  density  7Vtr 

2  x  1018 

cm-3 

SOA  internal  loss  ajnt 

5 

cm-1 

Injection  current  density  J 

3.4  x  103 

A/cm2 

Frequency  spacing  Q./  27r 

8.6 

GHz 

3)  Evaluate  the  field  vector  E(z  +  Ax)  by  solving  Eq.  (38) 
from  z  to  z  +  Ax  while  using  the  values  of  No  and  AiV/c 
obtained  in  steps  1  and  2,  according  to 

E(z  +  Ax)  =  exp  {  (1-fa)9WMl-a..AT 

exp  {T2m+i  [AAfc(x)]Ax}  E(z).  (46) 


A.  Model  Validation 


where 


-Psat(Ao) 


hu>oS 

r(N0)TgN(N0): 


(42) 


and  where  Ck  is  the  discrete  autocorrelation  function  of  E, 
namely 


M 

Ck=  J2  En+kEl  (43) 

nss-M 


in  which  we  assume  En  =  0  for  \n\  >  M.  The  expression  of  C 
in  the  case  M  =  2  is 
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where  we  used  C  /c  =  CJ[ ,  as  can  be  readily  verified  by  inspect¬ 
ing  Eq.  (43). 

The  numerical  integration  of  the  coupled-equations  involves 
a  three-step  procedure  for  the  transition  from  x  to  x  +  Ax,  given 
the  held  vector  E(z).  These  are: 

1)  Find  the  value  of  N0(z)  by  solving  Eq.  (39); 

2)  Extract  the  carrier  density  vector  A N(z)  as  in  Eq.  (41); 


In  this  section  we  test  the  accuracy  of  the  proposed  multi¬ 
wavelength  propagation  model  against  the  results  obtained  by 
integrating  the  full  model’  space-time  equations  (3)  and  ( 14).  To 
this  end  we  consider  a  QW  SOA,  characterized  by  the  following 
logarithmic  functional  dependence  of  the  gain  coefficient  on 
carrier  density  [6] 

g{N)  =  g0  log  ,  (47) 

where  go  is  a  gain  parameter  and  NtI  is  the  carrier  density 
required  for  transparency.4  The  expansion  of  the  gain  function 
is  in  this  case  g(N)  ~  g(N0)  +  g ^  (No) AN,  with 

g(N0)  =  g0  log  (J^)  >  9n(N0)  =  (48) 

The  physical  and  operational  parameters  of  the  SOA  are  listed 
in  Table  I  (we  note  that  the  SOA  is  operated  with  the  injection 
current  density  J  =  8.5 Jtr,  where  Jtr  =  ed(ANtI  +  BN t2r  + 
CJV/j )  is  the  injection  current  density  required  for  transparency). 
The  SOA  is  injected  with  a  three-wavelength  optical  signal  char¬ 
acterized  by  the  complex  envelope 

Ein  ( t )  =  ^fW[  e~mt  +  +  y/WTi  eMt  (49) 

with  W\  =  W- 1  =  —2  dBm,  and  Wo  =  —7  dBm.  For  this  set 
of  parameters  we  solved  the  coupled-mode  equations  (38),  (40), 
and  (41)  with  the  input  field  vector  E-in  =  [•  •  • ,  0,  \/W\ ,  yWo, 

4Of  course,  the  use  of  different  functional  forms  of  g(N),  for  instance  the 
more  accurate  three  parameter  expression  g(N)  =  go  In  [(TV  +  Ns)/(Ntr  + 
Ns )]  also  reported  in  [6],  is  fully  equivalent  in  terms  of  model  complexity. 
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Fig.  1 .  Intensity  (top  panel)  and  phase  (bottom  panel)  of  the  field  components 
Ek  versus  normalized  propagation  distance  z/L  for  the  displayed  values  of  k 
and  for  the  SOA  parameters’  values  in  Table  I.  Solid  curves  refer  to  the  coupled¬ 
mode  model,  while  circles  were  obtained  by  solving  the  space-time  equations 
of  the  full  model. 

y/W-i ,  0,  •  •  -]f  ■  We  used  M  =  6  and  checked  that  larger  values 
of  M  yield  indistinguishable  results.  We  then  integrated  the  full 
model’s  equations  (3)  and  (14)  with  the  procedure  described  in 
[19],  and  extracted  the  coefficients  Ek(z)  from  the  numerical 
solution  Enum(z,  t )  according  to 

Ek  (z)  ~  ~  r/n  Emm  (z,  tykatdt,  (50) 
Jto 

where  by  tg  we  denote  any  time  at  which  the  system  achieved 
its  stationary  state.  The  results  are  shown  in  Fig.  1.  In  the  top 
panel  we  plot  by  solid  curves  the  intensities  of  the  coefficients 
Ek(z)  versus  the  normalized  propagation  distance  z/L  for  val¬ 
ues  of  k  ranging  between  k  =  —4  and  k  =  4.  By  circles  we  plot 
the  results  obtained  with  the  full  model.  The  excellent  accu¬ 
racy  of  the  coupled-mode  model  is  self-evident.  Interestingly, 
the  figure  shows  that  the  coupled-mode  model  is  accurate  in 
describing  the  formation  of  FWM  components  that  eventually 
(at  the  SOA  output)  exceed  some  of  the  input  components.  The 
center  panel  of  the  same  figure  shows  the  corresponding  phases 
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Fig.  2.  Intensities  of  field  components  Ek  at  the  SOA  output  \Ek(L)\2,  as 
obtained  by  using  the  coupled-mode  model,  for  increasing  values  of  M .  Each 
curve  corresponds  to  a  specific  value  of  k  and  hence  it  originates  at  \k\  =  M. 
The  top  panel  refers  to  the  set  of  parameters  listed  in  Table  I  and  used  in  Fig.  1, 
whereas  in  the  bottom  panel  the  optical  confinement  factor  was  increased  from 
10%  to  20%. 


of  the  field  coefficients  Ek  (more  precisely  the  solid  curves  are 
the  plot  of  Phase  [Ek(z)\  +  kflz/vg,  where  the  second  term 
accounts  for  the  fact  that  the  coefficients  Ek  (z)  characterize  the 
field  envelope  in  the  time  reference  delayed  by  z/vg).  The  lower 
panel  provides  a  direct  validation  of  the  condition  j  A Nk  |  <C  No, 
which  was  assumed  in  the  perturbation  analysis.  The  top  circles 
show  the  time-averaged  carrier  density  (in  a  logarithmic  scale) 
extracted  from  the  full  model  results.  Although  this  quantity  is 
described  in  the  coupled-mode  model  by  Nq  +  AiVo,  we  com¬ 
pare  it  with  No,  which  is  shown  as  a  solid  curve.  The  excellent 
agreement  between  the  two  plotted  quantities  is  a  proof  of  the 
condition  ANg  <C  No.  Lower  curves  and  markers  are  plots  of 
\ANk  |  for  values  of  |fc|  ranging  between  1  and  4  (we  recall  that 
|AiVfcj  =  |  A N-k  |).  The  plot  shows  that  the  absolute  value  of  the 
coefficients  A Nk  is  more  than  two  orders  of  magnitude  smaller 
than  Nq,  thus  confirming  the  condition  A  A?/.  <C  No,  also  for 
|fc|  >  0.  Moreover,  the  excellent  agreement  between  the  full 
model  and  the  coupled-mode  model  validates  the  perturbation 
approach  for  arbitrary  magnitude  of  the  perturbations. 

We  stress  that  the  coupled-mode  model  offers  considerably 
greater  computational  efficiency,  as  compared  to  the  full  time 
domain  model.  In  the  specific  example  of  Fig.  1,  the  integration 
time  of  the  coupled-mode  equations  was  more  than  two  orders 
of  magnitude  smaller  than  the  integration  time  required  by  the 
space-time  equations,  using  in  both  cases  in-house  developed 
MATLAB  routines  run  on  the  same  workstation. 

In  Fig.  2  we  illustrate  the  dependence  of  the  coupled-mode 
model’s  results  on  the  number  of  field  coefficients  that  are 
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Fig.  3.  Dual-pumped  SOA-based  PSA’s  gain  versus  the  relative  phase  of 
the  input  signal  <j>8  introduced  in  Eq.  (55).  The  SOA  parameters  used  in  the 
numerical  computation  are  those  given  in  Table  I,  the  input  pump  powers  were 
settoVFp1  =  Wp2  =  —2  dBm,  and  the  input  signal  power  to  Ws  —  —22  dBm. 
The  solid  curve  refers  to  the  coupled-mode  model  with  M  =  4  (larger  values  of 
M  yield  indistinguishable  results),  while  the  circles  were  obtained  by  integrating 
the  space-time  equations  of  the  full  model.  The  dashed  curve  shows  the  results 
obtained  with  the  coupled-mode  model  by  propagating  only  the  pump  and  signal 
components,  namely  by  setting  M  =  1 . 

considered.  In  the  top  panel  we  plot  the  output  intensities 
|F4(L)|2  evaluated  by  solving  the  coupled-mode  equations  for 
increasing  values  of  M,  with  each  curve  corresponding  to  a 
different  value  of  k.  Since  the  accounting  for  the  frequency 
component  E *  dictates  that  M  >  |Zc|,  the  curve  referring  to 
Ej.  originates  at  M  =  \k\.  The  plot  shows  that  in  the  numer¬ 
ical  example  considered  here  the  results  of  the  coupled-mode 
equations  for  the  component  E *  become  accurate  (that  is,  the 
corresponding  curve  in  the  figure  becomes  flat)  for  M  exceeding 
\k\  by  a  one  or  two  units.  However,  it  should  be  pointed  out  that 
the  convergence  to  the  correct  result  is  affected  by  the  specific 
SOA  parameters’  value  and  may  be  slower.  This  is  shown  in 
the  lower  panel  of  the  Fig.  2,  where  the  same  curves  plotted  in 
the  top  panel  are  re-calculated  by  increasing  the  SOA  optical 
confinement  factor  from  10%  to  20%  and  by  leaving  the  other 
SOA  parameters  unchanged.  In  this  example,  it  can  be  seen  that 
using  M  <  6  may  yield  an  error  in  the  calculation  of  \Ei  (L)  |2 
up  to  a  factor  of  100. 

The  SOA  length  assumed  in  this  section  for  the  validation 
of  the  coupled-mode  model  is  intermediate  between  the  typical 
length  of  “ultra-short”  SOAs  ( L  <  0.5  mm)  and  that  of  “ultra- 
long”  SOAs  (L  >  1  mm).  We  note  that  while  the  accuracy  of  the 
coupled-mode  model  is  not  affected  in  the  case  of  shorter  SOAs, 
the  model  may  require  some  improvement  in  the  case  of  ultra- 
long  SOAs,  where  one  needs  to  take  into  account  the  spatial 
dependence  of  the  injection  current,  as  well  as  its  dependence 
on  the  local  carrier  density  [20]. 

IV.  Application  of  the  Model:  Dual-Pumped 
SOA-Based  PSA 

In  this  section  we  apply  the  coupled-mode  model  to  the  study 
of  the  dual-pumped  SOA-based  PSA  presented  in  [1 1] — [13]. 
The  goal  of  this  exercise  is  two-folded.  On  the  one  hand  we 
aim  to  show  that  the  phase-sensitive  gain  value  obtained  with 
the  coupled-mode  model  assuming  realistic  SOA  parameters  is 


consistent  with  the  experimentally  obtained  value.  On  the  other 
hand,  we  show  explicitly  that  by  restricting  the  coupled-mode 
model  to  the  pump  and  signal  components  only,  as  is  sometimes 
done  [21],  yields  significantly  incorrect  results  when  a  realistic 
dependence  of  the  amplifier  gain  on  carrier  density  is  used. 

The  waveform  at  the  input  of  a  PSA  of  the  kind  considered 
here  can  be  expressed  as 

Ein(t)  =  e^1  y/wie*1*  +  ei0°  V ,  +  e^-1  y/wl^eint, 

(51) 

where  by  W\  and  W- \  we  denote  the  optical  powers  of  the 
two  pumps  and  by  <j> \  and  <f>  \  their  absolute  phases.  The  field 
component  at  the  central  frequency  represents  the  input  sig¬ 
nal  component.  By  removing  in  all  components  the  immaterial 
average  phase  of  the  two  pumps 

0c  =  (0!  H-  0-i)/2,  (52) 

and  denoting  by 

<!>S  =  <A)  -  <I>C  (53) 

the  input  signal  phase  relative  to  0, ,  the  input  field  envelope  can 
be  expressed  as 

Ein(t)  =  ei^p  sJWle~im  +  ei<t>p  ^  , 

(54) 

where  <pp  =  (</>i  —  <^>_i  )/2.  We  further  note  that  the  effect  of 
is  limited  to  introducing  an  immaterial  time  shift  tp  = 
and  hence  it  can  be  safely  set  to  0V  =  0.  We  therefore  solve  the 
space-time  equations  using  the  following  input  waveform, 

Ein(t)  =  y/W^e-mt+ei*‘Vws  +  y/W^eint,  (55) 

and  the  coupled-mode  equations  with  the  input  field  vector 

Ein  =  [••■,  0,  y/W^,  €*•  sjws,  yf Wp, “  0,  ■■•]*•  (56) 

The  key  quantity  that  characterizes  the  performance  of  the  PSA 
under  scrutiny  is  the  dependence  of  the  signal  gain  on  the  phase 
(f)s.  In  Fig.  3  we  plot  the  gain  Gs((f>s)  =  \ES(L)\2 /Ws  (in  deci¬ 
bels)  as  a  function  of  <f)s ,  where  in  the  case  of  the  full  model  the 
term  Es  ( L )  is  extracted  from  the  numerical  solution  Enum  (L,  t) 
according  to  Eq.  (50)  with  z  =  L.  The  SOA  parameters  used  in 
the  numerical  example  are  those  given  in  Table  I.  The  input 
pump  powers  were  set  to  Wp1  =  Wp2  =  —2  dBm,  and  the  in¬ 
put  signal  power  to  the  much  smaller  value  Ws  =  —  22  dBm. 
The  solid  curve  is  obtained  by  integrating  the  coupled-mode 
model  with  M  =  4,  whereas  the  circles  refer  to  the  space- 
time  model.  The  excellent  agreement  between  the  coupled-mode 
model  and  the  space-time  model,  like  in  the  previous  section, 
is  self-evident.  The  thin  dashed  curve  in  Fig.  3  shows  the  result 
obtained  with  the  coupled-mode  model  by  including  only  the 
pump  and  signal  field  components,  that  is  by  using  M  —  1.  The 
plot  shows  that  the  neglect  of  high-order  FWM  products  yields 
higher  gain  values  and  a  lower  phase  dependent  gain. 

We  now  proceed  to  compare  the  results  of  the  coupled-mode 
model  with  the  experimental  data  reported  in  ref.  [13],  wherein 
all  the  details  concerning  the  experiment  can  be  found. 

The  comparison  between  the  theory  and  the  data  is  performed 
by  looking  at  the  dependence  of  the  PSA  gain  on  the  input  signal 
relative  phase  cj>s  of  the  kind  shown  in  Fig.  3.  To  this  end,  we 
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note  that  the  numerical  values  of  the  SOA  parameters  listed  in 
Table  I  are  those  estimated  for  one  of  the  devices  tested  in  [13], 
and  specifically  for  the  device  with  a  phase  dependent  gain  of 
about  6.3  dB.  Some  of  the  parameters  were  directly  taken  from 
the  device  geometry,  like  the  active  region  length  and  width. 
The  active  region  thickness  was  taken  to  be  65  nm  to  reflect 
the  fact  that  the  active  region  of  the  device  consisted  of  10 
InP/InGaAsP  QWs  of  6.5  nm  width  each.  Some  of  the  param¬ 
eters  were  estimated  by  independent  analysis,  like  for  instance 
the  value  of  T  ~  10%,  which  was  extracted  from  the  mode  pro¬ 
file  given  by  a  finite-difference-method-based  electromagnetic 
solver.  The  value  of  the  transparency  current  density  used  in 
simulations  is  Jth  =  0.43  kA/cm2,  versus  the  measured  value 
of  about  1  kA/cm2 .  This  implies  an  injection  efficiency  of  about 
43%,  which  includes  the  loss  caused  by  lateral  current  spread¬ 
ing  in  the  waveguide  region  and  the  carrier  trapping  efficiency 
in  the  QWs.  The  values  for  the  gain  coefficient  go  and  for  the 
coefficients  A,  B  and  C  are  those  typical  for  the  InP/InGaAsP 
MQW  active  region  of  the  SOA  used  in  the  experiment  [6], 
[13].  The  waveguide  internal  loss  was  also  set  to  the  typical 
value  of  dint  =  5  cnT1  measured  in  good  quality  devices.  The 
linewidth  enhancement  factor  a  =  5  is  also  within  the  typical 
range  of  values  for  devices  of  this  type.  With  these  parameters, 
the  model  predicts  about  59  dB  of  linear  gain,  whereas  the  mea¬ 
sured  linear  gain  of  the  device  under  test  was  about  48.5  dB 
[13].  The  10  dB  difference  can,  however,  be  safely  attributed  to 
gain  compression  induced  by  ASE — neglected  in  the  model  but 
significant  for  zero  input  in  devices  with  high  linear  gain — and 
to  thermal  effects  within  the  waveguide. 

The  comparison  between  the  theoretical  curve  shown  in  Fig.  3 
and  the  data  requires  some  care.  More  specifically,  since  the  way 
in  which  the  experiment  is  performed  does  not  allow  a  precise 
measurement  of  the  input  signal  power,  in  this  comparison  we 
look  at  the  output  signal  power  (rather  than  at  the  gain).  In 
the  logarithmic  scale,  this  simply  requires  adding  a  bias  to  the 
theoretical  PSA  gain.  The  relative  signal  phase  cf>s  defined  in 
Eq.  (53)  is  also  not  known  directly.  Indeed  in  the  experiment, 
the  phase  of  the  input  signal  is  controlled  by  passing  the  signal 
through  a  tuning  region  prior  to  the  injection  into  the  SOA.  The 
phase  change  of  the  signal  is  proportional  to  the  square  root 
of  the  current  /0  injected  in  the  tuning  region.  A  calibration 
procedure  allowed  to  establish,  for  between  1  and  2  mA, 
the  relation  tf>o  —  (f>\  =  7r X\/E  with  X  —  1  mA~l/2.  As  a  re¬ 
sult,  the  relation  between  the  signal  relative  phase  <ps  and  the 
measured  control  current  1^,  <j>s  =  ^X\fL  +  41 b ,  obtained  us¬ 
ing  the  definition  of  <ps  in  Eq.  (53),  contains  the  unknown  bias 
(f>b  =  (4> i  —  <^>_i)/2,  which  also  needs  to  be  added  to  the  theo¬ 
retical  signal  phase  in  order  to  facilitate  the  comparison  between 
the  theory  and  the  data. 

In  Fig.  4  we  plot  the  output  signal  power  versus  the  square- 
root  injection  current  yJT$.  The  markers  refer  to  the  experimen¬ 
tal  results,  and  they  were  taken  from  [13,  Fig.  26(a)].  The  solid 
curve  is  obtained  from  the  curve  plotted  in  Fig.  3  by  shifting  the 
vertical  and  the  horizontal  axes,  so  as  to  fit  the  data,  as  discussed 
in  the  previous  paragraph.  The  agreement  between  theory  and 
experiment  is  self  evident.  We  remark  that  such  good  agree¬ 
ment  cannot  be  obtained  by  using  the  conventional  three-wave 
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Fig.  4.  Measured  relationship  among  the  signal  power,  the  signal  phase  and 
the  square  root  of  the  current  applied  to  the  phase  tuner  (see  ref.  [13]),  together 
with  the  curve  of  Fig.  3,  suitably  re-centered  to  fit  the  data. 

model  (corresponding  to  the  case  M  =  1  in  our  theory  and  used 
in  ref.  [21])  where  the  high-order  harmonics  generated  by  the 
nonlinear  interaction  between  the  signal  and  the  two  pump  are 
neglected,  unless  unrealistic  values  for  the  device  parameters 
are  assumed. 

V.  Conclusion 

To  conclude,  we  derived  a  couple-mode  model  for  multi-wave 
mixing  in  SOAs  characterized  by  arbitrary  functional  depen¬ 
dencies  of  the  recombination  rate  and  material  gain  on  carrier 
density.  The  model  takes  into  account  the  frequency  dependence 
of  the  material  gain,  as  well  as  all  orders  of  the  waveguide  dis¬ 
persion,  and  accommodates  input  fields  consisting  of  arbitrary 
combinations  of  multiple  frequency  components.  We  showed 
that  the  conventional  approach  assuming  a  limited  number  of 
generated  FWM  components  gives  inaccurate  results  when  two 
waveforms  of  similar  intensities  are  injected  into  the  SOA.  In 
this  case,  our  model  gives  highly  accurate  results  if  a  sufficient 
number  of  generated  components  are  taken  into  account,  as  we 
showed  by  direct  comparison  with  full  time-domain  simulations. 
We  applied  the  coupled-mode  model  to  studying  the  operation 
of  a  recently  demonstrated  dual-pumped  PSA  based  on  an  inte¬ 
grated  QW  SOA  [1 1]— [13],  and  showed  that  the  outcome  of  the 
model  is  consistent  with  the  experimental  results. 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%  Three_input_waves_versus_z .m  %%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  This  file  produces  |E_k|A2  in  dB  versus  z/L,  and  Phase (E_k)  versus  z/L 
%  In  agreement  with  Fig.  1  of  the  paper 

clear;  clc;  close  all; 


%%%%%%%  Main  settings  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

WdBml  =  -2;  %  Input  power  at  frequency  w_0  -  Om  /  dBm 

WdBm2  =  -2;  %  Input  power  at  frequency  w_0  +  Om  /  dBm 

WdBmO  =  -7;  %  Input  power  at  frequency  w_0  /  dBm 

Ndz  =  1000;  %  Number  of  slices  in  the  device  modeling 

M  =  10;  %  Total  number  of  frequency  components  is  2M+1 

Mplot  =3;  %  The  plot  will  incude  2Mplot  +  1  frequency  components  around  w_0 .  N.B.  Mplot  <=  M 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%%%  INPUT  DATA  %%%%%%%%%% 

A  =  0*le6;  %  Linear  recombination  coefficient  /  sA-l 

B  =  0.3e-10;  %  Bimolecular  recombination  coefficient  /  cmA3/s 

C  =  3.3e-29;  %  Auger  recombination  coefficient  /  cmA3/s 

Gamma  =0.09;  %  Confinement  factor; 

lambdaO  =  1561e-9;  %  Wavelength  /  m 

D2  =  -8000;  %  Waveguide  dispersion  /  (fs/nm/cm), 

D2  =  D2*le-6;  %  Waveguide  dipersion  /  (s/m/ cm) 

beta2  =  -  1/ (2*pi*3e8 ) *lambda0A2*D2 ;  %  Waveguide  dipersion  /  (sA2/cm) 

alpha  =5;  %  Linewidth  enhancement  factor 

w  =  2e-4;  %  Width  of  the  waveguide  /  cm 

d  =  65e-7;  %  Active  region  tickness  /  cm 

L  =  0 . 1 ;  %  SOA  length  /  cm 

gNO  =  1800;  %  Gain  coefficent  /  cmA-l 

Ntr  =  2el8;  %  Transparency  carrier  density  /  cmA-3 

aint  =  0 . 5*le3*le-2 ;  %  Waveguide  internal  loss  /  cmA-l 

Ntrw  =  Ntr*exp (aint/ (Gamma*gN0) ) ;  %  Waveguide 

%  transparency  carrier  density  /  cmA-3 

e  =  1. 602e-19;  %  Electron  charge  /  C 

E  =  1.2825e-19;  %  Photon  energy  at  1.55  micrometers  /  J 

Jtr  =  (A*Ntrw+B*NtrwA2+C*NtrwA3 ) *e*d;  %  Waveguide  transparency  injection 

disp ([' Current  density  for  transparency  Jtr  =  num2str ( Jtr/1000) ,  'kA/cmA2  ']) 

%  current,  A/cmA2 

W1  =  10A (WdBml/10 ) *le-3;  %  Input  power  at  frequency  w_0  -  Om  /  W 

W2  =  10A (WdBm2/10 ) *le-3;  %  Input  power  at  frequency  w_0  +  Om  /  W 

W0  =  10A (WdBm0/10 ) *le-3;  %  Input  power  at  frequency  w_0  /  W 

lambdal  =  lambdaO  +  0.84e-9;  %  Wavelength  of  pump  2  (longer) 

lambda2  =  lambdaO  +  0.7e-9;  %  Wavelength  of  pump  1  (shorter) 
fl  =  3e8/lambdal; 
f2  =  3e8/lambda2; 

Om  =  2*pi* ( f2-f 1) /2 ;  %  2  Om  =  angular  frequency  shift  /  sA-l 
J  =  8 . 5* Jtr; %4500*0 . 8 ;  %  Injection  current  density  /  (A/cmA2) 
disp ([' Current  density  J  =  num2str ( J/1000 ) ,  '  kA/cmA2']) 

vg  =  le2*3e8/3 . 6; 

%%%%%%%%  END  OF  INPUT  DATA  %%%%%%%%%%%%%%% 

%%%%%%%%  Linear  gain  %%%%%%%%%%%%%%%%%%%%%%% 

Nin  =  fzero(@(N)  (A*N+B*NA2+C*NA3 ) - J/ (e*d) ,  Ntr);  %  Carrier  density  @  J  with  zero  optical  input 

Glin  =  exp ( (Gamma*gN0*log (Nin/Ntr)  -  aint)*L); 

disp ([' Linear  gain  Glin  =  num2str (10*logl0 (Glin) ) ,  '  dB ' ] ) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%  Preparing  to  loop  %%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
dz  =  L/Ndz;  %  Width  of  the  slice  /  cm 
Z  =  0:dz:L; 
count Z  =  1; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%  Field  initialization  % 

Nf  =  M;  %%%  Number  of  right/left  frequencies.  Total  is  2Nf+l 

N2pl  =  2*Nf +1 ; 

EE_in  =  zeros  (N2pl, 1) ; 

EE_in (Nf+1-1)  =  sqrt(Wl); 

EE_in (Nf+1  )  =  sqrt(WO); 

EE_in (Nf+1+1)  =  sqrt(W2);  %%%  Input  pump  2 

EE  =  EE_in; 

EEz  =  zeros (length (EE_in) , length ( Z) ) ;  EEz(:,l)  =  EE_in; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


for  z  =  Z  (2  : end) 
waitbar (z/L) 


%%%%%%  Working  point  %%%%%%%%%%%%%%%%% 

NO  =  fzero(@(x)  J/ (e*d)  - (A*x+B*xA2+C*xA3)  -  Gamma/ (E*w*d) *gN0*log (x/Ntr) * (EE ' *EE) 
RNO  =  A*N0+B*N0A2+C*N0A3; 

tau  =  1/ (A+2*B*N0+3*C*N0A2) ;  %  Spontaneus  lifetime  /  s 
Nin  =  NO; 

countZ  =  countZ  +1; 
gO  =  gN0*log (NO/Ntr) ; 

gN  =  gNO/NO;  %  Dif f erentiall  gain  at  reference  carrier  density  /  cmA2 
Psat  =  E* (w*d) / ( Gamma* gN* tau) ;  %  Local  saturation  power  /  W 
Pstim  =  E* (w*d) / (Gamma*gO  )*RN0; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


,  Nin) ;  %%%%  N_0  of  the  paper 


%%%  Evaluation  of  DN  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

EcorrE  =  conv (EE, conj (EE (end: -1 : 1) ) ,  ' same ' ) ; 

Cvec  =  EcorrE; 

Cvec(Nf+l)  =  0; 

Ck  =  ( [EcorrE (Nf+1 : -1 : 1 ); zeros (Nf, 1 )]) ; 

EcorrE  =  conv (EE, conj (EE (end: -1 : 1) ) ) ; 

Ck  =  EcorrE (2*Nf+l: -1:1) ; 

MatDNl  =  eye(2*Nf+l)  -  li*tau*Om*diag ( [Nf : -1 : -Nf ]  ); 

MatDN2  =  toeplitz (Ck) /Psat; 

DN  =  -  tau*RN0* (  MatDNl  +  MatDN2  ) A ( -1 ) *Cvec/Pstim; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%  Evolution  of  E  %%%%%%%%%%%%%%%% 

Emat  =  1/2* (1-li* alpha) *Gamma*gN* toeplitz ( [DN (Nf+1 : -1 : 1 ) ; zeros (Nf , 1 ) ] ) ; 

EE  =  expm (  (  ( ( l-li*alpha) * Gamma* gO-aint  ) /2*eye (2*Nf +1 )  +  Emat  ) *dz  )*EE; 
EEz  (:, countZ)  =  EE; 

DDN ( : , countZ)  =  DN; 

NOz  (countZ)  =  NO; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


end 


%%%%%%%  Compute  phase  of  Ek  %%%%%%%%%%%%%%%% 
for  nn  =  1:  (2*Nf+l) 
kk  =  nn  -  (Nf+1); 

phiz (nn, : )  =  angle (  EEz (nn, : ) . *exp (-li*kk*Om*Z/vg) 

end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


figure ( 1) ; 

plot (Z/L, 10*logl0 (abs (EEz (Nf+l-Mplot : Nf +l+Mplot , : ) ) . A2/le-3) , ' - ' , ' linewidth ' , 2) ; 
WmaxdBm  =  max (max (10*logl0 (abs (EEz (Nf+l-Mplot :Nf+l+Mplot, : ) ) . A2/le-3) ) ) ; 
ylim([-10,  WmaxdBm+1] ) ; 
xlabel ( ' z/L ' ) ;  ylabel ( ' I E_k | A2  (dBm) ' ) 

figure (2) ; 

plot (Z/L, phiz (Nf+l-Mplot : Nf +l+Mplot , : ) , ' - ' , ' linewidth ' , 2 ) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%  PSA_gain_versus_phi_s .m  %%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  This  file  produces  signal  gain  in  dB  versus  signal  relative  phase  phi_s 
clear;  clc;  close  all; 

%%%%%%%%%%  Selection  of  input  and  output  %%%%%%%%% 

Min  =0;  %  Input  frequency:  Min  *  Om 
Mout  =  2;  %  Output  frequency:  Mout  *  Om 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%  Main  settings  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

WdBml  =  -2-3;  %  Input  pump  power  at  frequency  w_0  -  Om  /  dBm 
WdBm2  =  -2;  %  Input  pump  power  at  frequency  w_0  +  Om  /  dBm 

WdBmO  =  min ( [WdBml, WdBm2 ] )  -  20;  %  Input  signal  power  at  frequency  w_0  /  dBm 

Ndz  =  100;  %  Number  of  slices  in  the  device  modeling 

M  =  5;  %  Total  number  of  frequency  components  is  2M+1 

phi_s_vec  =  (0: 0. 001: 1) *pi;  %  Signal  phase  values  phi_s 

%%%%  Computation  time  increases  by  increasing  Ndz 

%%%%  and  by  -3increasing  the  resolution  of  phi_s_vec 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%  INPUT  DATA  %%%%%%%% 

A  =  0*le6;  %  Linear  recombination  coefficient  /  sA-l 

B  =  0.3e-10;  %  Bimolecular  recombination  coefficient  /  cmA3/s 

C  =  3.3e-29;  %  Auger  recombination  coefficient  /  cmA3/s 

Gamma  =  0.098;  %  Confinement  factor; 

lambdaO  =  1561e-9;  %  Wavelength  /  m 

D2  =  -8000;  %  Waveguide  dispersion  /  (fs/nm/cm), 

D2  =  D2*le-6;  %  Waveguide  dipersion  /  (s/m/ cm) 

beta2  =  -  1/ (2*pi*3e8 ) *lambda0A2*D2 ;  %  Waveguide  dipersion  /  (sA2/cm) 

alpha  =5;  %  Linewidth  enhancement  factor 

w  =  2e-4;  %  Width  of  the  waveguide  /  cm 

d  =  65e-7;  %  Active  region  tickness  /  cm 

L  =  0.1;  %  SOA  length  /  cm 

gNO  =  1800;  %  Gain  coefficent  /  cmA-l 

Ntr  =  2el8;  %  Transparency  carrier  density  /  cmA-3 

aint  =  0 . 5*le3*le-2 ;  %  Waveguide  internal  loss  /  cmA-l 

Ntrw  =  Ntr*exp (aint/ (Gamma*gN0) ) ;  %  Waveguide 

%  transparency  carrier  density  /  cmA-3 

e  =  1.602e-19;  %  Electron  charge  /  C 

E  =  1.2825e-19;  %  Photon  energy  at  1.55  micrometers  /  J 

Jtr  =  (A*Ntrw+B*NtrwA2+C*NtrwA3) *e*d;  %  Waveguide  transparency  injection 

disp ([' Current  density  for  transparency  Jtr  =  ',  num2str ( Jtr/1000) ,  'kA/cmA2  ']) 

%  current,  A/cmA2 

W1  =  10 A (WdBml/10 ) *le-3;  %  Input  power  at  frequency  w_0  -  Om  /  W 

W2  =  10 A (WdBm2/10 ) *le-3;  %  Input  power  at  frequency  w_0  +  Om  /  W 

W0  =  10 A (WdBm0/10 ) *le-3;  %  Input  power  at  frequency  w_0  /  W 

lambdal  =  lambdaO  +  0.84e-9;  %  Wavelength  of  pump  2  (longer) 

lambda2  =  lambdaO  +  0.7e-9;  %  Wavelength  of  pump  1  (shorter) 
fl  =  3e8/lambdal; 
f2  =  3e8/lambda2; 

Om  =  2*pi* ( f2-f 1) /2 ;  %  2  Om  =  angular  frequency  shift  /  sA-l 
J  =  8 . 5* Jtr ; %4500*0 . 8 ;  %  Injection  current  density  /  (A/cmA2) 
disp ([' Current  density  J  =  ',  num2str ( J/1000 ) ,  '  kA/cmA2 ' ] ) 

vg  =  le2*3e8/3 . 6; 

%%%%%%%%  END  OF  INPUT  DATA  %%%%%%%%%%%%%%% 

%%%%%%%%  Linear  gain  %%%%%%%%%%%%%%%%%%%%%%% 

Nin  =  fzero(@(N)  (A*N+B*NA2+C*NA3 ) - J/ (e*d) ,  Ntr);  %  Carrier  density  @  J  with  zero  optical  input 

Glin  =  exp ( (Gamma*gN0*log (Nin/Ntr)  -  aint)*L); 

disp ([' Linear  gain  Glin  =  ',  num2str (10*logl0 (Glin) ) ,  '  dB ' ] ) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%  Preparing  to  loop  %%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
dz  =  L/Ndz;  %  Width  of  the  slice  /  cm 
Z  =  0 : dz : L; 
count_phi  =  1 ; 

Nseed  =  Nin; 

Pt  =  zeros (1, Ndz) ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%%  Field  initialization  % 

Nf  =  M;  %%%  Number  of  right/left  frequencies.  Total  is  2Nf+l 


N2pl  =  2*Nf+l; 

EE_in  =  zeros (N2pl, 1) ; 
EE_in (Nf+1-1)  =  sqrt(Wl); 
EE_in (Nf+1+1)  =  sqrt (W2 ) ; 
EE  -  EE_in; 


EEz  =  zeros (length (EE_in) , length ( Z) 
E_gain  =  sqrt (WO); 
nphi  =  length (phi_s_vec) ; 
gs  =  zeros (1, nphi) ; 
phi_vec  =  zeros (1 , nphi) ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%  Input  pump  2 
EEz  (  : , 1)  =  EE_in; 


for  phi  -  phi_s_vec 
waitbar (phi/pi) 


%%%%%%  Per-phase  value  initialization  %%%%%%%%%%%%%%%%% 

EE  -  EE_in; 

EE (Nf  +  1  +  Min)  =  sqrt (WO) *exp (-li*phi) ;  %%%  Input  signal 

count Z  =  1; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


z  =  Z (2 : end) 

%%%%  Working  point  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
if  count_phi  ==  1  %%%  Solve  the  fzero  equation 

Nin  =  fzero (@ (x)  - (A*x+B*xA2+C*xA3 ) + J/ (e*d)  -  Gamma*gNO*log (x/Ntr) * (EE ' *EE ) / (E*w*d)  ,  Nseed)  . 

Nseed  =  Nin; 

Pt(countZ)  =  EE ' *EE ; 
count Z  =  countZ  +1; 
else  %%%  Use  the  the  look-up  table 

[XXX,  Pt_index]  =  min (abs (Pt_vec-EE ' *EE) ) ; 

Nin  =  Nin_vec (Pt_index) ; 

end 

NO  =  Nin; 

RNO  =  A*N0+B*N0  A2+C*N0 A  3; 

tau  =  1/ (A+2*B*N0+3*C*N0A2) ;  %  Spontaneus  lifetime  /  s 
Nin  =  NO; 

gO  =  gN0*log (NO/Ntr) ; 

gN  =  gNO/NO;  %  Dif f erentiall  gain  at  reference  carrier  density  /  cmA2 
Psat  =  E* (w*d) / ( Gamma* gN*tau) ;  %  Local  saturation  power  /  W 
Pstim  =  E* (w*d) / (Gamma*gO  )*RN0; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%  Evaluation  of  DN  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

EcorrE  =  conv (EE, conj (EE (end: -1 : 1 ) ) , ' same ' ) ; 

Cvec  =  EcorrE; 

Cvec(Nf+l)  =  0; 

EcorrE  =  conv (EE, conj (EE (end: -1 : 1 ) ) ) ; 

Ck  =  EcorrE (2*Nf+l: -1:1) ; 

MatDNl  =  eye(2*Nf+l)  -  li*tau*Om*diag (Nf : -1 : -Nf ) ; 

MatDN2  =  toeplitz (Ck) /Psat; 

DN  =  -  tau* RNO* (  MatDNl  +  MatDN2  ) A ( -1 ) *Cvec/Pstim; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%%%  Evolution  of  E  %%%%%%%%%%%%%%%% 

Emat  =  1/2* (1-li* alpha) *Gamma*gN* toeplitz ( [DN (Nf+1 : -1 : 1 ) ; zeros (Nf , 1 ) ] ) ; 

EE  =  expm (  (  ( ( l-li*alpha) *Gamma*gO-aint  ) /2*eye (2*Nf +1 )  +  Emat  ) *dz  )*EE; 
if  count_phi  —  1 

E_gain  =  exp ( ( ( l-li*alpha) *Gamma*gO-aint) /2*dz) *E_gain; 

end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


end 

%%%  Build  the  the  look-up  table  %%%%%%%%%%%%%%%%%% 
if  count_phi  —  1 
NPt  =  le4 ; 

Nseed  =  Ntr; 

Pt_vec  =  zeros ( 1, NPt+1) ; 

Nin_vec  =  zeros (1 , NPt+1 ) ; 
countPt  =  1; 

for  Pt  =  (0 :NPt) /NPt*2*max (Pt) 

Pt_vec (countPt)  =  Pt; 

Nin_vec (countPt)  =  fzero (@(x)  - (A*x+B*xA2+C*xA3) +J/ (e*d)  -  Gamma*gN0*log (x/Ntr ) *Pt/ (E*w*d)  , 

Nseed  =  Nin_vec (countPt) ; 
countPt  =  countPt  +  1; 

end 


Nseed) 


end 


phi_vec (count_phi)  =  phi; 

Ws  =  abs (EE (Nf+l+Mout) A2) ; 
gs (count_phi)  =  10*logl0 (Ws/WO) ; 
count_phi  =  count_phi  +  1; 

end 

Gsat  =  10*logl0 (abs (E_gain) A2/W0) ; 

disp ([' Saturated  gain  Gsat  =  ',  num2str (Gsat) ,  '  dB']) 

figure ( 3) ;  plot (phi_vec, gs, '  k ' , [0 , pi] ,  Gsat* [1,1], ' — g', ' linewidth ' , 2 ) 
xlim ( [0  pi] ) ; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%  PSA_NF_estimate.m  %%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  This  file  produces  noise  figure  in  dB  versus  signal  relative  phase  phi_s 
clear;  clc;  close  all;  %rng ( 1  shuffle ') ; 

%%%%%%%%%%  Selection  of  input  and  output  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
Min  =0;  %  Input  frequency:  Min  *  Om 
Mout  =0;  %  Output  frequency:  Mout  *  Om 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%  Main  settings  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
WdBml  =  -2;  %  Input  pump  power  at  frequency  w_0  -  Om  /  dBm 

WdBm2  =  -2;  %  Input  pump  power  at  frequency  w_0  +  Om  /  dBm 

WdBmO  =  -24  +  min ( [WdBml , WdBm2 ]) ;  %  Input  signal  power  at  frequency  w_0  /  dBm 

Ndz  =  100;  %  Number  of  slices  in  the  device  modeling 

M  =  5;  %  Total  number  of  frequency  components  is  2M+1 

phi_s_vec  =  (0 . 005 : 0 . 01 : 1 . 005) *pi;  %  Signal  phase  values  phi_s 

%%%%  Computation  time  increases  by  increasing  Ndz 

%%%%  and  by  increasing  the  resolution  of  phi_s_vec 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%  INPUT  DATA  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

A  =  0*le6;  %  Linear  recombination  coefficient  /  sA-l 

B  =  0.3e-10;  %  Bimolecular  recombination  coefficient  /  cmA3/s 

C  =  3.3e-29;  %  Auger  recombination  coefficient  /  cmA3/s 

Gamma  =0.098;  %  Confinement  factor; 

lambdaO  =  1561e-9;  %  Wavelength  /  m 

D2  =  -8000;  %  Waveguide  dispersion  /  (fs/nm/cm) , 

D2  =  D2*le-6;  %  Waveguide  dipersion  /  (s/m/ cm) 

beta2  =  -  1/ (2*pi*3e8 ) *lambda0A2*D2 ;  %  Waveguide  dipersion  /  (sA2/cm) 

alpha  =5;  %  Linewidth  enhancement  factor 

w  =  2e-4;  %  Width  of  the  waveguide  /  cm 

d  =  65e-7;  %  Active  region  tickness  /  cm 

L  =  0.1;  %  SOA  length  /  cm 

gNO  =  1800;  %  Gain  coefficent  /  cmA-l 

Ntr  =  2el8;  %  Transparency  carrier  density  /  cmA-3 

aint  =  0*0 . 5*le3*le-2 ;  %  Waveguide  internal  loss  /  cmA-l 

Ntrw  =  Nt r* exp (aint/ (Gamma* gNO) ) ;  %  Waveguide 

%  transparency  carrier  density  /  cmA-3 

e  =  1.602e-19;  %  Electron  charge  /  C 

E  =  1.2825e-19;  %  Photon  energy  at  1.55  micrometers  /  J 

Jtr  =  (A*Ntrw+B*NtrwA2+C*NtrwA3) *e*d;  %  Waveguide  transparency  injection 
disp ([' Current  density  for  transparency  Jtr  =  ',  ... 

num2str ( Jtr/1000) ,  'kA/cmA2  ']) 

%  current,  A/cmA2 

W1  =  10A (WdBml/10) *le-3;  %  Input  power  at  frequency  w_0  -  Om  /  W 

W2  =  10A (WdBm2/10) *le-3;  %  Input  power  at  frequency  w_0  +  Om  /  W 

W0  =  10A (WdBm0/10 ) *le-3;  %  Input  power  at  frequency  w_0  /  W 

lambdal  =  lambdaO  +  0.84e-9;  %  Wavelength  of  pump  2  (longer) 

lambda2  =  lambdaO  +  0.7e-9;  %  Wavelength  of  pump  1  (shorter) 
fl  =  3e 8 /lambdal ; 
f 2  =  3e8/lambda2 ; 

Om  =  2*pi* ( f2-f 1) /2 ;  %  2  Om  =  angular  frequency  shift  /  sA-l 
J  =  8 . 5* Jtr ; %4500*0 . 8 ;  %  Injection  current  density  /  (A/cmA2) 
disp ([' Current  density  J  =  num2str ( J/1000 ) ,  '  kA/cmA2 ' ] ) 

vg  =  le2*3e8/3 . 6; 

Bw  -  le3; 

%%%%%%%%  END  OF  INPUT  DATA  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%  Linear  gain  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
Nin  =  fzero(@(N)  (A*N+B*NA2+C*NA3) -J/ (e*d) ,  Ntr);  %  Carrier  density  @  J 
%with  zero  optical  input 

Glin  =  exp ( (Gamma*gN0*log (Nin/Ntr)  -  aint)*L); 
disp ([' Linear  gain  Glin  =  num2str (10*logl0 (Glin) ) ,  '  dB ' ] ) 
disp ([' Unsaturated  carrier  density  =  num2str (Nin) ,  '  mA-3']) 
disp ([' Transparency  carrier  density  =  num2str (Ntr ) ,  '  mA-3']) 

nspin  =  Nin/ (Nin-Ntr) ; 

%  nspin  =  (Nin/ (Nin-Ntr )) A . 2 ; 
disp (['Input  nsp  =  num2str (nspin) ] ) 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%  Preparing  to  loop  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
dz  =  L/Ndz;  %  Width  of  the  slice  /  cm 
Z  =  0 : dz : L; 
count_phi  -  1 ; 

Nseed  =  Nin; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%%  Field  initialization  % 

Nf  =  M;  %%%  Number  of  right/left  frequencies.  Total  is  2Nf+l 

N2pl  =  2*Nf+l; 

EE_in  =  zeros  (N2pl, 1) ; 

EE_in (Nf+1-1)  =  sqrt(Wl); 

EE_in (Nf+1+1)  =  sqrt(W2);  %%%  Input  pump  2 

nphi  =  length (phi_s_vec) ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

Es  =  zeros ( 1, nphi) ; 
phi_vec  =  zeros (1 , nphi) ; 

Pt  =  zeros ( 1, Ndz) ; 

%%%%%%  Per-phase  value  initialization  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
EE  =  EE_in; 

EE(Nf+l+Min)  =  sqrt(WO);  %%%  Input  signal 
E_gain  =  sqrt(WO); 
count Z  =  1; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%Part  1:  BUILDING  A  LOOKUP  TABLE  FOR  N  vs.  THE  AVERAGE  INTENSITY%%%%%%%% 


for  z  =  Z (2 : end) 

%%%%  Working  point  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%  Solve  the  fzero  equation 

Nin  =  fzero (@ (x)  - (A*x+B*xA2+C*xA3) + J/ (e*d)  ... 

-  Gamma*gNO*log (x/Ntr) * (EE ' *EE) / (E*w*d)  ,  Nseed); 

Nseed  =  Nin; 

Pt(countZ)  =  EE'*EE; 
count Z  =  countZ  +1; 

NO  =  Nin; 

RNO  =  A*N0+B*N0 A2+C*N0A3; 

tau  =  1/ (A+2*B*N0+3*C*N0A2) ;  %  Spontaneus  lifetime  /  s 
Nin  =  NO; 

gO  =  gN0*log (NO/Ntr) ; 

gN  =  gNO/NO;  %  Dif f erentiall  gain  at  reference  carrier  density  /  cmA2 
Psat  =  E* (w*d) / (Gamma*gN*tau) ;  %  Local  saturation  power  /  W 
Pstim  =  E* (w*d) / (Gamma*gO  )*RN0; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%  Evaluation  of  DN  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

EcorrE  =  conv (EE, conj (EE (end: -1 : 1) ) ,  ' same ' ) ; 

Cvec  =  EcorrE; 

Cvec(Nf+l)  =  0; 

EcorrE  =  conv (EE, conj (EE (end: -1 : 1 ) ) ) ; 

Ck  =  EcorrE (2*Nf+l: -1:1) ; 

MatDNl  =  eye(2*Nf+l)  -  li*tau*Om*diag (Nf : -1 : -Nf ) ; 

MatDN2  =  toeplitz (Ck) /Psat; 

DN  =  -  tau* RNO* (  MatDNl  +  MatDN2  )A(-l)*Cvec/Pstim; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%  Evolution  of  E  %%%%%%%%%%%%%%%% 

Emat  =  1/2* (l-li*alpha) *Gamma*gN*toeplitz ( [DN (Nf+1 :-l :1) ; zeros (Nf ,1) ] ) ; 
EE  =  expm ( ( ( ( l-li*alpha) * Gamma* gO-aint  ) /2*eye (2*Nf+l) +  Emat  ) *dz  ) *EE; 
E_gain  =  exp (  ( (l-li*alpha) *Gamma*gO-aint  ) *dz/2) *E_gain; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


end 

NPt  =  le4 ; 

Nseed  =  Ntr; 

Pt_vec  =  zeros ( 1, NPt+1) ; 

Nin_vec  =  zeros (1 , NPt+1 ) ; 
countPt  =  1; 

for  Pt  =  (0 :NPt) /NPt*2*max (Pt) 

Pt_vec (countPt)  =  Pt; 

Nin_vec (countPt)  =  fzero (@(x)  - (A*x+B*xA2+C*xA3) +J/ (e*d)  ... 

-  Gamma*gN0*log (x/Ntr) *Pt/ (E*w*d)  ,  Nseed); 

Nseed  =  Nin_vec (countPt) ; 
countPt  =  countPt  +  1; 

end 

Gsat  =  abs (E_gain) A2/W0 ; 

GsatdB  =  10*logl0 (Gsat) ; 

disp ([' Saturated  gain  Gsat  =  num2str (GsatdB) ,  '  dB']) 

%%%%Part  2:  EVALUATING  THE  NOISLESS  OUTPUT%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

for  phi  -  phi_s_vec 


waitbar (phi/pi) 


%%%%%%  Per-phase  value  initialization  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

EE  =  EE_in; 

EE(Nf+l+Min)  =  sqrt (WO) *exp (-li*phi) ;  %%%  Input  signal 

count Z  =  1; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

for  z  =  Z (2 : end) 

%%%%  Working  point  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%  Use  the  the  look-up  table 

[XXX,  Pt_index]  =  min (abs (Pt_vec-EE ' *EE) ) ; 

Nin  =  Nin_vec (Pt_index) ; 

NO  =  Nin; 

RNO  =  A*N0+B*N0 A2+C*N0A3; 

tau  =  1/ (A+2*B*N0+3*C*N0A2) ;  %  Spontaneus  lifetime  /  s 
Nin  =  NO; 

gO  =  gN0*log (NO/Ntr) ; 

gN  =  gNO/NO;  %  Dif f erentiall  gain  at  reference  carrier  density/cmA2 
Psat  =  E* (w*d) / ( Gamma* gN* tau) ;  %  Local  saturation  power  /  W 
Pstim  =  E* (w*d) / (Gamma* gO  )*RN0; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%  Evaluation  of  DN  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

EcorrE  =  conv (EE, conj (EE (end: -1 : 1 ) ) , ' same ' ) ; 

Cvec  =  EcorrE; 

Cvec(Nf+l)  =  0; 

EcorrE  =  conv (EE, conj (EE (end: -1 : 1) ) ) ; 

Ck  =  EcorrE (2*Nf+l: -1:1) ; 

MatDNl  =  eye(2*Nf+l)  -  li*tau*Om*diag (Nf : -1 : -Nf ) ; 

MatDN2  =  toeplitz (Ck) /Psat; 

DN  =  -  tau* RNO * (  MatDNl  +  MatDN2  ) A ( -1 ) *Cvec/Pstim; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%  Evolution  of  E  %%%%%%%%%%%%%%%% 

Emat  =  1/2* (l-li*alpha) *Gamma*gN*toeplitz ( [DN (Nf+1 : -1 : 1 ) ; zeros (Nf,l) ] ) ; 

EE  =  expm (  (  ( ( l-li*alpha) *Gamma*gO-aint  ) /2*eye (2*Nf +1 )  +  Emat  ) *dz  ) *EE 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

phi_vec (count_phi)  =  phi; 

Es (count_phi)  =  EE (Nf +l+Mout) ; 
count_phi  =  count_phi  +  1; 

end 

gs  =  10*logl0 (abs (Es) . A2/W0 ) ; 

figure (3);  plot (phi_vec, gs,  ' k ' ,  [ 0 , pi ] ,  GsatdB* [1, 1] ,  ' — c ' , ' linewidth ' , 2 ) 

nsrout  =  zeros ( 1, nphi ); %  Initialization  of  the  output  noise  to  signal  ratio 
nsroutl  =  zeros (1 , nphi) ; 

nsrin  =  (E*Bw/2)/W0;  %  Input  noise  to  signal  ratio 

%%%%%%Part  3:  MONTE  CARLO  RUNS  FOR  FOR  THE  SNR%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


for  Nloop  =  1:5000 

%%%%%%%%  Linear  gain  %%%%%%%%%%%%%%%%%%%%%%% 

Nin  =  fzero(@(N)  (A*N+B*NA2+C*NA3) -J/ (e*d) ,  Ntr) ;  %  Carrier  density 
%  @  J  with  zero  optical  input 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%  Preparing  to  loop  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

count_phi  =  1 ; 

Nseed  =  Nin; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

Els  =  zeros (1 , nphi) ; 
for  phi  -  phi_s_vec 
waitbar (phi/pi) 

%%%%%%  Per-phase  value  initialization  %%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
EE  -  EE_in; 

EE(Nf+l+Min)  =  sqrt (WO) *exp (-li*phi) ;  %%%  Input  signal 

EE  =  EE  +  sqrt (Bw*E/2) * (randn (N2pl, 1) +li*randn (N2pl,  1) ) ;  %  input 

%  field  with  zero  point  noise  over  bandwidth  Bw  (power  spectral 

%  density  of  half  a  photon  per  quadrature) 

count Z  =  1; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

for  z  =  Z  (2  : end) 

[XXX,  Pt_index]  =  min (abs (Pt_vec-EE ' *EE) ) ; 

Nin  =  Nin_vec (Pt_index) ; 

NO  =  Nin; 

RNO  =  A*N0+B*N0A2+C*N0A3; 

tau  =  1/ (A+2*B*N0+3*C*N0A2) ;  %  Spontaneus  lifetime  /  s 
%  nsp  =  NO/ (NO-Ntr) ; 

%  nsp  =  (NO/ (NO-Ntr) ) A . 2; 
nsp  =  1; 

Nin  =  NO; 

gO  =  gN0*log (NO/Ntr) ; 

gN  =  gNO/NO;  %Dif ferential  gain  at  reference  carr.  density/cmA2 
Psat  =  E* (w*d) / (Gamma*gN*tau) ;  %  Local  saturation  power  /  W 
Pstim  =  E* (w*d) / (Gamma*gO ) *RN0 ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%  Evaluation  of  DN  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

EcorrE  =  conv (EE, conj (EE (end: -1:1) ) , ' same ' ) ; 

Cvec  =  EcorrE; 

Cvec(Nf+l)  =  0; 

EcorrE  =  conv (EE, conj (EE (end:-l :1) ) ) ; 

Ck  =  EcorrE (2*Nf+l:-l :1) ; 

MatDNl  =  eye(2*Nf+l)  -  li*tau*Om*diag (Nf : -1 : -Nf ) ; 

MatDN2  =  toeplitz (Ck) /Psat; 

DN  =  -  tau* RNO* (  MatDNl  +  MatDN2  ) A ( -1 ) *Cvec/Pstim; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%%%  Evolution  of  E  %%%%%%%%%%%%%%%% 

Emat  =  1/2* (l-li*alpha) *Gamma*gN . . . 

* toeplitz ( [DN (Nf+1 : -1 : 1 ) ; zeros (Nf , 1 ) ] ) ; 

EE  =  expm ( ( ( ( l-li*alpha) *Gamma*gO-aint)  ... 

/2*eye (2*Nf+l) +  Emat) *dz  ) *EE  ... 

+  sqrt (E*Bw* ( (2*nsp-l ) *Gamma*gO+aint) *dz/2)  ... 

* (randn (N2pl, 1) +li*randn (N2pl, 1) ) ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


end 

phi_vec (count_phi)  =  phi; 

Els (count_phi)  =  EE (Nf+l+Mout) ; 
count_phi  =  count_phi  +  1; 

end 

nsrout  =  nsrout* (Nloop-1) /Nloop  ... 

+  real (  (Es-Els)  . *exp (-li*angle (Es) ) )  . A2  ... 

. / real (Es . *exp (-li* angle (Es) ) ) . A2 /Nloop; 

%  Updated  output  noise  to  signal  ratio 
nsroutl  =  nsroutl* (Nloop-1) /Nloop  . . . 

+  imag ( (Es-Els) . *exp (-li*angle (Es) )). A2  ... 

. /real (Es . *exp (-li*angle (Es) ) ) . A2/Nloop; 

%  Updated  output  ortogonal  noise  to  signal  ratio 

figure (4) ; 

plot (phi_vec, 10*logl0 (nsrout/ nsrin) , ' o — r ' ,  ... 

phi_vec, 10*logl0 (nsroutl /nsrin)  ,  ' o — b ' ,  ... 

[0,  pi],  [3,  3],'—  c',  [0,  pi],  [0,  0],'— g'); 

xlim( [0,pi] ) ; 

disp(['End  of  iteration  #  ',  num2str (Nloop) ] ) 

end 

nsrin  =  (E*Bw/2)/W0;  %  Input  noise  to  signal  ratio 
figure (5) ; 

plot (phi_vec, 10*logl0 ( (nsrout+nsroutl ) / (2*nsrin) )  ,  ' o — k' , 
[0,  pi],  [3,  3],'—  b\  [0,  pi],  [0,  0],'— g'); 
xlim ( [0,pi] ) ; 


